Objective Metasurfaces enable flexible manipulation of electromagnetic waves with planar and integrated architectures, showing great potential for compact optical imaging and sensing systems. Nevertheless, full-wave electromagnetic simulation becomes a critical bottleneck for large-aperture high-numerical-aperture (NA) metalenses. Conventional finite-difference time-domain (FDTD) simulation demands prohibitive memory and runtime resources as device dimensions increase. Many reported acceleration techniques rely on periodic hypothesis or symmetry constraints, which cannot faithfully capture complex inter-element near-field interactions. Existing acceleration approaches including local periodic approximation (LPA) field-stitching suffer severe accuracy degradation under high-NA conditions because strong aperiodic near-field coupling between adjacent meta-units cannot be properly captured. Widely adopted domain-decomposition schemes still lack quantitative trade-off rules among accuracy, peak memory and computation time for two-dimensional large-aperture metalenses. Practical applications such as astronomical imaging, head-mounted virtual-reality displays and large-aperture optical sensing systems require millimeter-scale metalens devices, yet reliable full-wave performance assessment remains difficult under conventional general-purpose computing infrastructure. This work aims to resolve the memory bottleneck of FDTD full-wave simulation for millimeter-scale high-NA metalenses, establish quantitative performance trade-offs for sub-regional simulation, and provide a practical simulation workflow for high-fidelity performance evaluation of large-area metasurface devices.

Methods Quadratic scaling relationships between FDTD computational overhead and metalens physical area are quantitatively derived under uniform mesh discretization which is widely adopted in commercial time-domain electromagnetic solvers. The error boundaries of the LPA field-stitching technique are systematically assessed by comparing stitched near-field distributions against full-wave reference results across varying NA values. Normalized global error is adopted as the core figure-of-merit to quantify near-field discrepancies by integrating vector-field residual energy over the complete metalens output aperture. A sub-regional simulation framework based on perfectly-matched-layer (PML) absorbing boundaries is proposed. Each sub-domain consists of an inner effective region and surrounding transition-buffer units to suppress artificial boundary-truncation perturbations and reconstruct realistic aperiodic electromagnetic environments. Heat-map-based parameter sweeping is performed to quantify how effective-unit number and transition-buffer-unit number influence sub-domain global error, single-subdomain runtime and memory consumption. Resource extrapolation is conducted for a 1-mm-diameter metalens to predict peak memory and serial total runtime. The proposed sub-regional scheme is validated against full-wave simulation and LPA stitching on a 25-μm-aperture high-NA metalens. Near-field vector electric-field components and far-field focusing characteristics are compared for quantitative evaluation.

Results and Discussions FDTD grid count, runtime memory and simulation time scale quadratically with metalens diameter under uniform meshing. A millimeter-scale metalens requires approximately 7.4 TB of peak memory for direct full-wave simulation, which exceeds the capacity of ordinary computing workstations and most small-scale research clusters. LPA stitching delivers acceptable accuracy only for metalenses with slowly-varying meta-unit geometries. Global error rises rapidly with increasing NA; considerable amplitude distortion and phase deviation emerge over the whole aperture at NA = 0.95. Sharp geometric variations of meta-units near the lens periphery excite prominent longitudinal electric-field components and inter-element scattering, which further aggravate field reconstruction errors and makes LPA unsuitable for high-precision evaluation. Introducing adequate transition-buffer units significantly mitigates boundary-induced disturbances inside the effective sub-domain region. Increasing unit numbers reduces simulation error but raises per-subdomain memory and time cost. Optimized sub-domain parameters yield 4.8% sub-domain global error while compressing peak memory requirement from 7.4 TB down to 1.9 GB for the 1-mm-diameter device. Serial total runtime may increase owing to accumulated computation over a huge number of independent sub-domains. However, sub-domains require no mutual data communication, which allows straightforward distributed parallel deployment across multiple workstations. Parallel computation can greatly compress total wall-clock time without raising peak-memory requirements. Validation on the 25-μm-aperture metalens shows that sub-regional simulation reconstructs vector near-field components and far-field focal patterns highly consistent with full-wave results, achieving 1.13 % global error. By contrast, LPA stitching fails to reproduce transverse and longitudinal vector field components and produces large overall deviation.

Conclusions The quadratic resource-scaling law reveals intrinsic computational burdens for large-aperture metalens full-wave simulation. LPA exhibits clear application limits caused by ignored aperiodic near-field coupling at high NA. The PML-based sub-regional simulation framework converts TB-level global simulation tasks into GB-level sub-domain computations. Peak memory is drastically reduced at the cost of moderately increased serial runtime. Quantitative trade-off rules among accuracy, memory and runtime guide the selection of effective-region and transition-buffer sizes. Different from symmetry-restricted decomposition methods, this sub-regional strategy imposes no rotational-symmetry constraints and supports general metasurface devices including beam-splitters and holographic elements. Distributed parallel capability further improves practical feasibility. This framework offers an efficient engineering tool for performance assessment, and can serve as verification support for large-aperture metasurface inverse-design workflows. Serial computation overhead remains the main limitation; future improvements may incorporate symmetry boundary conditions, GPU acceleration and artificial-intelligence surrogate models to shorten total computation time.