Efficient $C^1$ Bernstein Quasi-Trefftz Discretization for Heterogeneous High-Frequency Helmholtz Problems

• Bernstein–Bézier coordinates support exact geometric smoothness constraints and let heterogeneous material coefficients enter through local algebra rather than a special-function basis.Differentiation, edge restriction, normal differentiation, and degree elevation are sparse coefficient operations.Source: §1 · p. 1
• On full-row-rank two-triangle patches, the method reduces p2 + p + 1 conforming coordinates to 2p + 1 equation-adapted coordinates before global coupling.The global system therefore acts on trace-sized local vectors rather than full polynomial B-nets.Source: §1 · pp. 1–2
• The construction separates geometry-only C1 reduction, coefficient- and operator-dependent residual reduction, and global coupling through traces, cell residuals, and boundary residuals.The residual reduction is formed inside the conforming space, before neighboring patches interact.Source: §1 · pp. 1–2
• The method extends quasi-Trefftz and embedded-Trefftz algebraic reductions by applying them after strong C1 compatibility and expressing them directly in Bernstein coordinates.This representation is also used for heterogeneous coefficients, nonzero forcing, local compression, and reduced global assembly.Source: §1 · p. 2
• The implementation combines explicit C1 synthesis, QR-based local reduction, projected coefficients, graph-residual and physical-energy estimates, and batched reduced-coordinate assembly.For the exact-circle scattering test with 200 macro-patches and 3400 active unknowns, optimized preprocessing takes less than one second.Source: §1 · p. 2
• The computational study examines conformity, conditioning, heterogeneous-medium convergence, resolution thresholds, coefficient-projection sensitivity, penetrable scattering, and an Airy turning-point transition.The tests include resolved and polluted high-frequency regimes at fixed κh/p.Source: §1 · p. 2

2. Bernstein–Bézier construction of the local C1 space

• B-coefficients are constrained to form the classical spline space, and the differential equation is used only after that exact C1 space has been formed.The B-net is the collection of Bernstein coefficients over the domain points.Source: §2 · p. 3

2.1. Smoothness conditions across one interior edge

• Across an interior edge, matching the edge coefficients gives continuity, while a first-order relation between adjacent coefficient rows supplies the additional C1 conditions.For an isolated two-triangle join, the construction imposes p + 1 value conditions and p first-order smoothness conditions.Source: §2.1 · p. 3

2.2. Derivative and trace maps in B-coordinates

• Directional differentiation is sparse on neighboring B-coefficients and lowers polynomial degree; traces use the edge row, while normal derivatives also use the first adjacent row.The same normal is used on both adjacent triangles to orient the derivative trace consistently.Source: §2.2 · p. 4

2.3. Assembly of the macro-patch spline space

• The assembled smoothness matrix depends only on geometry and degree, while multi-triangle dimensions depend on smoothness and residual ranks rather than a universal two-triangle formula.The construction forms the conforming space first and imposes PDE residual equations in those coordinates.Source: §2.3 · p. 5

3.1. What quasi-Trefftz means in the present method

• When the degree-(p-2) residual map has full row rank, a two-triangle C1 patch retains 2p+1 local coordinates instead of quadratic-growth conforming polynomial coordinates.The low-order interior residual is eliminated before patches communicate; the global unknowns are the surviving local coordinates.Source: §3.1 · pp. 5–6

3.2. Residual moments by integration by parts

• Residual moments are assembled by integration by parts, requiring point values of A and η but never derivatives of A.Triangle-edge boundary contributions are retained; they cancel for a continuous exact medium and a C1 patch field.Source: §3.2 · p. 6

3.3. Adaptive algebraic-degree projection of nonpolynomial coefficients

• Adaptive algebraic degree projection selects local polynomial approximations of A and η, while the global graph residual uses the original coefficients and source.The experiments use εA = 10^-9; tightening it changes reported errors below displayed digits.Source: §3.3 · p. 7

3.4. Coefficient approximation and residual-map perturbation

• The local kernel is stable when coefficient-induced residual-map perturbations remain below the smallest retained singular value, with rotation controlled relative to that gap.A complete perturbation norm includes the conormal boundary contribution, not only an L2(K) coefficient indicator.Source: §3.4 · pp. 8–9

3.5. Affine quasi-Trefftz patch space for nonzero forcing

• In heterogeneous tests, every patch retained exactly 2p+1 homogeneous coordinates and projected forcing moments matched to roundoff.For nonzero forcing, the source changes only the affine lift; globally coupled unknowns remain homogeneous coordinates.Source: §3.5 · p. 9

3.6. The surviving quasi-Trefftz consistency defect

• The oversampled true residual splits into coefficient/source approximation error and Trefftz truncation error, with the latter containing only modes above degree p-2.If projected medium and source are exact and the polynomial space contains a local exact solution, both terms vanish and coupling becomes trace-only.Source: §3.6 · p. 10

4. Intrinsic graph-residual coupling

• The graph residual combines scaled value/conormal mismatches with a cell term measuring the consistency defect that survives local reduction.On full-row-rank two-triangle patches, local coordinates fall from O(p^2) to 2p + 1; an exact Trefftz basis makes the cell term vanish.Source: §4 · p. 11

4.1. Residual optimality and discrete injectivity

• If the assembled residual map is injective on the homogeneous reduced space, the least-squares problem has a unique solution and a Hermitian positive-definite normal matrix.Source: §4.1 · p. 11
• The discrete solution is an exact best approximation in the graph norm, while physical-energy reliability follows by transferring residual and Cauchy-mismatch control.The transfer explicitly removes locally annihilated degree-(p −2) residual modes before global assembly.Source: §4.1 · p. 12, §4.2 · p. 12

4.2. Physical-energy reliability at high frequency

• The low-order AAD defect is controlled by coefficient/source approximation, while wave resolution is governed separately by κh/p.For exact-coefficient quasi-Trefftz errors, the low-order defect vanishes.Source: §4.2 · p. 13
• The physical-energy transfer separates coefficient/source error from the unresolved high-order residual tail, whose bound gains a factor h_P/(p+1).For exact polynomial media with the physical residual represented in the degree-p channel, both additional terms vanish.Source: §4.2 · pp. 13–14
• The reliability bound assumes continuous impedance stability and Cauchy lifting; it does not establish wavenumber-uniform impedance constants for general heterogeneous anisotropic media.Its discrete resolution factor is χ_h = max_P κh_P/p; the continuous constants may still depend on frequency, domain, and coefficient geometry.Source: §4.2 · pp. 14–15
• The coefficient-approximation contribution has a separately controllable indicator, but the exact-coefficient residual matrix needed for it is unavailable for nonpolynomial operators.The practical operator-relevant stopping test is suggested; the reported experiments use the coefficient-space test.Source: §4.2 · p. 15

5.1. Explicit C1 synthesis on a two-triangle patch

• Explicit Bernstein continuation fixes the common edge and first adjacent row, giving an explicit coordinate system for the two-triangle C1 spline space.The constructed space agrees with a rank-revealing nullspace computation to principal-angle accuracy near roundoff.Source: §5.1 · p. 15

5.2. Local quasi-Tre

• Adjoint QR extracts the local residual kernel; on full-row-rank two-triangle patches its dimension is 2p + 1.Direct elimination has condition numbers growing from about 10^4 at p = 6 to about 10^12 at p = 14, while the reduced basis condition number remains about 1.4–1.6.Source: §5.2 · p. 16
• Reference-triangle maps and affine-local geometry factors support the same Bernstein–Bézier smoothness construction on arbitrary noncongruent triangles.Coefficient projection and residual-moment assembly are batched using precomputed reference projection systems.Source: §5.2 · p. 16

5.4. Reduced global graph assembly

• Global traces, conormal fluxes, and the volume residual are assembled after contraction with the reduced basis, leaving only neighboring patches coupled through interface rows.The resulting normal matrix follows the macro-patch adjacency graph; the volume channel tests the complement of residual directions already annihilated locally.Source: §5.4 · p. 16

5.5. Equivalence and timing checks

• At p = 12–16, explicit C1 synthesis is about 100–170 times faster than rank-revealing nullspace formation; batched global assembly is about 3.8 times faster at p = 8–10.Restricting the volume channel to the surviving residual complement cuts volume work by factors of 3.5–5.8; optimized stages agree with unreduced predecessors to roundoff-level matrix differences.Source: §5.5 · p. 17

6. Numerical experiments

• Inside each two-triangle patch, the polynomial pieces join strongly in C1; neighboring patches couple through value and conormal-flux mismatches in the graph residual.The experiments use unstructured unit-square triangulations paired into nonoverlapping two-triangle patches.Source: §6 · p. 17, §6.1 · p. 17, Figure 1 · p. 18

6.1. Strong C1 conformity on noncongruent patches

• Both the full spline basis and reduced basis satisfy the smoothness conditions to roundoff, with dimensions p^2 + p + 1 and 2p + 1, respectively, under full row rank.Table 1 reports worst mismatches over the mesh and the range of nonzero-singular-value condition numbers of the scaled local smoothness matrix.Source: §6.1 · p. 17, Table 1 · p. 18

6.2. Conditioning of the reduced global system

• In the tested unstructured cases, the residual operator remains moderately conditioned, and diagonal scaling reduces the normal-matrix condition in every displayed case.Table 2 reports conditioning of the graph-residual system on unstructured macro-patch meshes.Source: §6.2 · p. 19, Table 2 · p. 19

6.3. A matrix-valued heterogeneous medium and unstructured refinement

• For κ=40 and p=8, convergence persists across noncongruent patches with varying triangle shapes, while each patch retains 17 active coordinates and the C1 construction remains unchanged.The test uses a nonpolynomial matrix-valued heterogeneous medium with spatially varying anisotropy.Source: §6.3 · p. 19, Figure 3 · p. 20

6.4. Fixed-mesh degree resolution and high wavenumber

• The fixed-mesh degree sweep shows rapid error decrease after κh/p reaches the local resolution threshold, with κ=40 reaching 7.04 × 10^-8 at p=15.At κ=80, the error drops sharply near κh/p≈1; κ=120 remains underresolved on this mesh through p=15.Source: §6.4 · pp. 19, 21, Figure 4 · p. 21
• At κ=80, relative L2 error falls from 3.65 (p=10) to 1.38 × 10^-3 (p=16), while κ=120 remains unresolved through p=15.A finer 256-patch run at κ=120, p=16 reaches relative L2 error 5.79 × 10^-2 with 8448 active coordinates.Source: Figure 4 · p. 21, §6.4 · p. 21

6.5. Pollution at fixed resolution

• Intermediate targets show fitted exponents near zero at κh/p≈0.6 and positive at ≈0.7, placing 0.7 in the transition between resolved and near-threshold regimes.The exponents at 0.6 are −0.03 and −0.18 for p=8 and 10; at 0.7 they are 0.28 and 0.65.Source: §6.5 · p. 21
• At κh/p≈0.5, endpoint error ratios from κ=20 to 110 are 0.74 (p=8) and 0.84 (p=10), versus 13.4 and 36.4 at κh/p≈0.8.The near-threshold log–log fits have descriptive exponents 0.95 and 1.78, respectively.Source: §6.5 · p. 21, Figure 5 · p. 22

6.6. Penetrable scattering with an exact NtD truncation

• At κ=40 and p=8, the exact-disk scattering run has a 4.77 s runtime and 1.04 × 10^-3 relative NtD boundary residual, with 17 coordinates per patch.The exact circular boundary avoids geometric boundary approximation error; the computed field exhibits refraction, focusing, and asymmetric interference.Source: §6.6 · p. 22, Figure 6 · p. 23

6.7. Sensitivity to the local coefficient projection

• Once coefficient approximation falls below discretization error, tightening the local projection tolerance leaves the field unchanged to displayed digits; the 8.13 × 10^-4 saturation is spatial discretization error.The operator-level stopping test is more selective when the retained singular gap is small and is intended for rank-sensitive regimes.Source: §6.7 · pp. 22, 24

6.8. An Airy turning-point transition beyond the uniformly positive regime

• The Airy test lies outside the uniformly positive-index assumptions: the exact field is evanescent on one side of the turning point and oscillatory on the other.For the sign-changing coefficient, the boundary least-squares term uses a fixed positive Robin weight.Source: §6.8 · p. 24
• The same conformity-first construction crosses the turning line cleanly, retaining trace-sized ranks 13, 17, and 21 for p=6, 8, and 10.The numerical study reports that the local quasi-Trefftz mechanism remains effective as the solution changes character across the domain.Source: §6.8 · p. 24, Table 6 · p. 24

Relation to Trefftz methods

• On full-row-rank two-triangle patches, removing p(p − 1) residual directions leaves 2p + 1 coordinates, and fixed-mesh accuracy improves rapidly once the local resolution threshold is crossed.The original strongly C1 space has p^2 + p + 1 coordinates; the reduced global dimension grows only linearly with p.Source: §7 · p. 24
• Unlike Taylor quasi-Trefftz methods, which annihilate a residual jet near a point [10], this method imposes projected residual moments across the macro-patch; unlike embedded Trefftz methods [15] [21], reduction follows C1 geometric restriction.This ordering is the source of the conformity-first structure.Source: §7 · p. 26

Computational implications

• Like plane-wave DG and related Trefftz-DG methods, global coupling uses traces, but this method retains the consistency defect in a graph functional and uses polynomials for heterogeneous coefficients and forcing.It avoids selecting propagation directions or replacing the local basis with special functions.Source: §7 · p. 26
• Local construction temporarily uses a quadratic-size conforming space, while sparse global coupling uses linear-size reduced coordinates and the macro-patch adjacency graph.Explicit C1 continuation, batching, and QR reduction keep preprocessing cost small in the tested range; geometry-dependent quantities can be cached across frequency sweeps on a fixed mesh.Source: §7 · p. 26
• Physical-energy stability constants remain problem dependent, with impedance and Cauchy-lifting constants depending on wavenumber, anisotropy, and domain; the analysis does not claim a wavenumber-uniform heterogeneous bound.The estimate separates these continuous stability factors from the discrete κh/p dependence.Source: §7 · p. 26
• The implementation uses two-triangle patches; larger aggregates may have topology-dependent smoothness dependencies or residual-rank loss, motivating further study of sparse continuation and patch-shape effects.The general multi-triangle algorithm assembles the smoothness matrix before forming the residual map in conforming coordinates.Source: §7 · p. 26

8. Conclusions

• Projected coefficients construct the local reduced space, while the original medium enters the global residual, separating coefficient approximation from discretization error.The perturbation analysis quantifies kernel rotation through the retained singular gap and separates represented high-order residual, low-order coefficient/source defect, and unresolved residual tail.Source: §8 · p. 26
• The method maintains strong C1 conformity to roundoff, achieves high-order convergence after the local resolution threshold, and remains stable in Airy turning-point and exact-circle scattering tests.Fixed-resolution experiments distinguish little frequency growth from a near-threshold polluted regime; the Airy test lies outside the positive-index assumptions used in the physical-energy estimate.Source: §8 · p. 26
• Dense work is concentrated in local preprocessing while the global system uses reduced coordinates; preprocessing is no longer the dominant cost in the representative scattering experiment.The authors identify repeated solves, parameter sweeps, and adaptive extensions as settings for which this separation is attractive.Source: §8 · p. 27

Data and code availability

• The source archive includes the unstructured mesh generator, local reduction and graph-assembly routines, scattering and Airy drivers, degree-sweep driver, figure scripts, and numerical data.Source: Data and code availability · p. 27

Declaration of generative AI and AI-assisted technologies in the writing process

• The author reports using ChatGPT for code prototyping and testing, numerical experiments, visualization, mathematical checking, and drafting and editing, and states that the manuscript was reviewed and verified.Source: Declaration of generative AI and… · p. 27