Comb Tensor Network Finite-Difference Time-Domain Algorithms for 2D and 3D Electromagnetic Simulations

1 Introduction

• Standard FDTD costs O(NtN^d) computation and O(N^d) memory, making fine discretizations and dense PML or material profiles expensive.Fine grids and dense spatial representations can suppress numerical dispersion and resolve geometry details, but raise memory and CPU demands.Source: §1 · p. 1
• Tensor-train and quantized tensor-train FDTD can reduce complexity on structured meshes and smooth solutions, but compression depends on tensor-site ordering and complex fields can compress poorly.Existing tensor-compressed approaches may also require storing or manipulating intermediate 2D or 3D arrays.Source: §1 · p. 2
• CTN tensorization aligns with multidimensional grid operators; its 3D structure differs from conventional QTT, while its 2D form is equivalent to QTT with seq(BF) ordering.The paper motivates CTN for directional coupling in Maxwell equations and anisotropic operators and boundary effects in FDTD with PML.Source: §1 · p. 2
• CTN-FDTD represents fields, permittivity, and auxiliary PML variables as comb tensor networks, applying finite-difference operators without forming full spatial arrays during time stepping.Local and global smoothing suppress rank inflation from repeated derivative updates and noise accumulation.Source: §1 · p. 2
• 35 trillion spatial grid points are tractable for CTN-FDTD in 3D, whereas standard FDTD would require 281 TB to store one 3D array.The paper reports reduced memory and CPU complexity while preserving accuracy and stability across 2D and 3D benchmarks.Source: §1 · p. 2

2 The Standard FDTD Algorithms with PML

• The formulation considers Maxwell equations in a lossless, unbounded domain with constant permeability, while allowing permittivity to vary across the domain.The computational grid includes the PML region and uses d = 2 or 3 spatial dimensions.Source: §2 · pp. 2–3

2.1 2D FDTD with split-field PML

• The 2D split-field PML formulation splits Ez into Ezx and Ezy and updates each N × N field quantity at staggered time steps.The described finite differences approximate first-order spatial derivatives; outside the PML, specified coefficient matrices take constant values.Source: §2.1 · p. 3

2.2 3D FDTD with ADE-PML

• The 3D ADE-PML formulation adds 12 auxiliary field quantities alongside six electric and magnetic field quantities, discretized as 18 N × N × N tensors.The update uses conductivity and shift parameters in PML regions and a discretized curl operator.Source: §2.2 · pp. 3–4
• 31 mode-3 tensors are required at minimum by the described 3D ADE-PML algorithm, including fields, auxiliary quantities, coefficients, and inverse relative permittivity.At N = 2048, the passage reports at least 2.13 TB of memory, and notes that this requirement becomes prohibitively high for large N.Source: §2.2 · pp. 3–4

3 CTN-FDTD Algorithms with PML

• CTN augmentation gives storage O(R^2 log N) and computation O(R^3N log N), where R is the CTN rank.The algorithms described cover 2D and 3D FDTD with PML, assuming N = 2^L.Source: §3 · p. 4
• The method assumes N = 2^L for both 2D and 3D simulations.The section presents operations mostly for 3D and briefly notes their adaptation to 2D.Source: §3 · p. 4
• The stated complexity is parameterized by CTN rank R, so its scaling depends on the rank as well as grid size.The passage gives distinct rank powers for storage and computation: R^2 and R^3, respectively.Source: §3 · p. 4

3.1 CTN representation of FDTD quantities

• A 3D CTN organizes tensor cores into a spine and branches, with the comb representation formed by joining the specified cores.The CTN is also described as closely related to extended TT or QTT-Tucker decompositions.Source: §3.1 · pp. 4–5
• The bond indices γ1 and γ2 connect the tensor cores and have dimensions R1 and R2, respectively.The indices i, j, and k denote grid points along the three spatial dimensions.Source: §3.1 · p. 5
• The 2D CTN is essentially equivalent to QTT using the specialized seq(BF) site ordering, whereas the 3D CTN is structurally different from QTT.In seq(BF), the two dimension index groups use opposite least-significant-bit to most-significant-bit orderings.Source: §3.1 · p. 5

3.2.1 Hadamard product

• Exact Hadamard products multiply CTN ranks; SVD-based truncation, optionally using zip-up on each TT branch, reduces them.The product is elementwise and applies to CTNs of the same dimensions.Source: §3.2.1 · p. 6

3.2.2 Summation

• CTN summation adds ranks, after which SVD-based truncation reduces them; scalar-weighted sums work similarly.Source: §3.2.2 · p. 7

3.2.3 Kronecker product of QTTs

• A Kronecker product of tensorized vectors uses their QTT cores as CTN branch cores and sets the spine cores to 1.The vectors are tensorized along specified spatial dimensions.Source: §3.2.3 · p. 7

3.2.4 Shift

• Finite-difference derivatives use shifts along a spatial mode; applying the QTT shift operator updates only the corresponding CTN branch cores.The shift maps each vector entry f_i to f_{i+1} and has an efficient QTT representation.Source: §3.2.4 · p. 7

• The 3D auxiliary arrays and excitation and field quantities are constructed as CTNs without forming full N × N × N tensors.The same principle applies to N × N tensors in 2D.Source: §3.3 · p. 7
• Separable auxiliary quantities have CTN representations requiring at most O(N) computation and peak storage; the example rank is at most 2Npml.The stated bound applies to a PML auxiliary quantity and extends to the other listed auxiliary quantities.Source: §3.3 · p. 7
• The considered inverse permittivity is a sum of a few separable components, giving O(N) construction costs; TT-cross-type algorithms can approximate more complex profiles.Source: §3.3 · p. 7
• Localized Gaussian, all-one, or standard-basis excitation vectors have CTN representations constructible with O(log N) computation and peak storage.Source: §3.3 · p. 8

3.4 Smoothing Techniques

• Updates can increase field ranks through accumulated numerical errors and complex fields in inhomogeneous permittivity, motivating two smoothing strategies.Source: §3.4 · p. 8

3.4.1 Local smoothing

• Local smoothing is applied after updating each field quantity and uses scalar products, sums, and shifts implemented on CTNs.Source: §3.4.1 · p. 8
• Local smoothing adds artificial dissipation through an even-order spatial derivative; this paper uses orders 2 and 6.The coefficient satisfies 0 ≤ η ≪ dt/ds^m and controls dissipation strength while maintaining update stability.Source: §3.4.1 · p. 8

3.4.2 Global smoothing

• Global smoothing also applies to inverse permittivity and uses forward DFT, low-pass filtering, inverse DFT, then SVD-based rank truncation.The filter may be Gaussian or Sigmoid; the DFTs and filter act branchwise without modifying the CTN spine cores.Source: §3.4.2 · pp. 8–9

• Initialization costs are O(dN) for time-independent auxiliary quantities and O(d log N) for time-dependent excitation; TT-cross can reduce the former if needed.Source: §3.5 · p. 9
• Each time step costs O(dR^3 log N) computation, dominated by SVD-based rank truncation, while storage is O(dR^2 log N).R is the maximum CTN rank across fields and time steps, and is not necessarily constant with N.Source: §3.5 · p. 9
• The CTN update and smoothing costs contrast with standard FDTD’s O(Nd) computation per time step and O(Nd) storage.Table 1 summarizes CPU and memory complexities for Nd grid points and Nt = O(N) time steps.Source: §3.5 · p. 10, Table 1 · p. 10
• Reconstructing fields for visualization adds O(N2) computation and memory, but the analysis excludes this optional postprocessing.In 3D, only N × N slices are reconstructed; in 2D, full N × N tensors are reconstructed.Source: §3.5 · p. 10

4 Numerical Results

• The examples demonstrate CTN-FDTD accuracy, efficiency, and scalability across 2D and 3D cases.Both CTN-FDTD and standard FDTD were implemented in Python.Source: §4 · p. 10
• All experiments ran on one CPU core of a Perlmutter node with two 64-core AMD EPYC 7763 processors.The node had 512 GB of 2133MHz DDR4 memory.Source: §4 · p. 10
• The Perlmutter node was located at NERSC in Berkeley.The experiments used one CPU node.Source: §4 · p. 10

4.1.1 Line Current Source

• At L = 15, standard FDTD requires 20× more CPU time and 105× more memory than CTN-FDTD for the line-current test.The comparison covers L = 8 to 15 at PPW=20; CTN-FDTD follows O(N log N) CPU and O(log N) memory scaling.Source: Figure 3 · p. 11, §4.1.1 · p. 11

4.1.2 Gaussian Pulse Initial Field

• With the Gaussian initial field, CTN-FDTD uses substantially less memory but beats standard FDTD in runtime only at L = 15.Observed maximum ranks for L = 9 to 15 are 63, 132, 159, 149, 172, 254, 343, and 247; the larger ranks yield scalings above the estimates.Source: Figure 4 · p. 12, §4.1.2 · p. 12

4.1.3 Point Current Source

• For the point-source test at fixed frequency, CTN-FDTD shows O(N log N) CPU and O(log N) memory scaling, beating standard FDTD in memory for L > 8 and CPU for L > 13.With PPW fixed at 7, its CPU time is not comparable to standard FDTD in the high-frequency regime.Source: §4.1.3 · p. 13
• For the point source, smoothing reduces ranks from near the theoretical maximum of 2L = 512 to below 70, while distorting field patterns to some extent.Sixth-order local smoothing gives the best observed tradeoff between accuracy and low-rankness.Source: §4.1.3 · p. 13, Figure 5 · p. 14

4.1.4 Photonic Crystal Waveguide

• In the photonic-crystal waveguide, ranks grow slightly but remain well below 2^10 = 1024, consistent with field confinement inside the waveguide.The test uses Gaussian global smoothing for 1/ϵr and second-order local smoothing for the fields.Source: §4.1.4 · p. 14

4.2.1 Plane-wave Excitation without Scatterer

• In homogeneous plane-wave propagation, CTN ranks remain small constants over time, while η = 0.0001 can produce larger ranks for theoretically zero Ex and Hx fields.The test uses f = 16 MHz, L = 10, and 6th-order local smoothing; Figure 10 plots field patterns, ranks, and memory.Source: §4.2.1 · p. 14, Figure 10 · p. 18
• At L = 15, CTN-FDTD handles 35 trillion grid points with projected 10 orders of magnitude less memory and 6 orders of magnitude faster computation than standard FDTD.Standard FDTD could only run up to L = 10, approximately 1 billion points; smoothing keeps CTN ranks constant as L increases.Source: §4.2.1 · pp. 14–15

4.2.2 Plane-wave Excitation with Scatterer

• The scatterer test uses a 200 m cubic inclusion with ϵr = 2, illuminated by two plane waves in a 1 × 1 × 1 km3 domain.The grid has L = 11 and the simulation runs 3500 time steps.Source: §4.2.2 · pp. 15, 17
• With a cubic scatterer, capping all CTN ranks at 64 completes 3500 time steps without visible field differences and captures complex field patterns.Global smoothing of 1/ϵr eliminates diffraction and reflection due to the sharp boundary; auxiliary-field ranks grow as near-zero amplitudes accumulate low-rank approximation error.Source: §4.2.2 · p. 17

• CTN-augmented 2D and 3D FDTD algorithms represent fields and permittivity with CTNs and use local and global smoothing to reduce ranks.Source: §5 · p. 18
• Numerical results demonstrate O(N log N) CPU and O(log N) memory complexities for d-dimensional domains with O(Nd) grid points across varied excitations, frequencies, and scatterers.Source: §5 · p. 18
• Future work includes parallelization, further reducing computational complexity, and improving stability for more complex scatterers and media.Source: §5 · p. 18