Accurate prediction of tensorial spectra using equivariant graph neural network

Domain: chemistry
Published: 2026-03-02
Editorial outcome: major_revision

Abstract
Optical spectroscopies provide a powerful means to probe light–matter interactions and complex electronic features that are crucial for the development and optimization of optoelectronic devices, where performance is closely tied to the underlying electronic spectrum. However, realistic modeling of tensor optical responses in materials remains computationally demanding and challenging. Here we introduce the Tensorial Spectra Equivariant Neural Network (TSENN), an equivariant graph neural network architecture that maps crystal structures directly to their full photon-frequency-dependent optical tensors. By encoding isotropic sequential scalar components and anisotropic sequential tensor components into spherical tensor representations, TSENN ensures symmetry-aware predictions consistent with crystalline symmetry constraints. Trained on frequency-dependent dielectric tensors of 1,432 bulk semiconductors, the model achieves a mean absolute error of 0.127, demonstrating its potential for efficient and general modeling of optical properties. Our framework opens new avenues for data-driven design of anisotropic optical responses to accelerate materials discovery for optoelectronic applications.

Introduction
Optical spectroscopies, which involve excitation of a material with external electromagnetic fields, provide a powerful tool for probing properties of materials. The two key frequency-dependent response functions involved here are the dielectric function and optical conductivity, which together form the foundation for interpreting a wide range of optical phenomena in solid-state systems [1–3]. Dielectric function captures the energy storage and dissipation mechanisms in the material, while optical conductivity characterizes the current response to varying external fields. Refractive index, extinction coefficient, and absorption coefficient also play important roles in specific applications such as the LEDs, integrated circuits, and solar cells [4–8]. Optical responses are often divided into various frequency regions, where the spectra carry fingerprints of distinct types of quasiparticles in materials [9–11].

Materials with low-symmetry crystal structures exhibit pronounced anisotropies which require tensors with non-zero off-diagonal components and a basis for new strategies for controlling their optical properties. For example, direction-dependent optical responses of van der Waals heterostructures can be tailored for polarization-sensitive photodetectors and developing thermal management schemes, among other applications [12, 13]. These heterostructures are also promising materials platforms for developing ultrafast lasers and artificial synaptic devices, and other transformational emerging optoelectronic technologies [14–16].

Graph neural networks (GNNs) have recently shown strong capabilities in capturing structure–property relationships and efficiently predicting complex functional properties of materials [17–20]. In particular, Euclidean neural networks (E(3)NN) [21–23], which explicitly incorporate spatial symmetry constraints into their architecture can, in principle, capture the full symmetry of a crystal. This enables high-fidelity predictions with relatively small training datasets (∼ 10^3 examples). E(3)NNs have demonstrated success in a variety of applications, including the prediction of interatomic potentials [24] and the construction of machine-learning Hamiltonians [25, 26].

Equivariant GNNs have also been extended to predict sequential material properties such as phonon density of states, optical conductivities, and absorption spectra [18, 27–29]. For instance, GNNOpt employs ensemble embeddings of aggregated node features to predict the trace of dielectric functions and derived properties, enabling the discovery of promising solar–cell materials and the estimation of quantum weights in materials [30]. Similarly, Optimate and its extension were initially trained on 9,915 semiconductors and insulators using dielectric functions within the independent-particle approximation (IPA), and subsequently fine-tuned on random-phase approximation (RPA) spectra to enhance their correspondence with experimental observations [31, 32]. In parallel, several works have proposed equivariant architectures for predicting tensor quantities such as dielectric, piezoelectric, and elastic tensors [33, 34]. For example, AnisoNet leverages spherical harmonic decomposition in the ℓ = 0 and ℓ = 2 channels to directly predict the dielectric tensor, facilitating high-anisotropy materials discovery through database screening [35], while DTNet adopts a transfer-learning approach that incorporates interatomic potentials as optimized latent embeddings to construct the 3×3 static dielectric constant using node features and atomic positions [36]. Notably, most existing studies focus exclusively on the dielectric constant. In contrast, the joint prediction of tensor structure and spectral behavior–specifically, the full frequency-dependent dielectric tensor–remains largely unexplored. Prior work has investigated tensorial spectra indirectly through equivariant machine–learning Hamiltonians, but the Brillouin-zone integrations involved continue to pose a computational bottleneck. Our approach addresses this gap by leveraging the expressive power of higher-order tensor channels within the equivariant GNNs to directly predict the full, frequency-dependent dielectric tensor from crystal structure, thereby providing a comprehensive framework for modeling anisotropic optical responses in crystalline materials.

Here we propose the Tensorial Spectra Equivariant Neural Network (TSENN) whose architecture is based on E(3)NN, which is equivariant under translation, rotation and inversion symmetries. These symmetry-aware capabilities of E(3)NN allow TSENN to model complex directional dependencies of functional properties of materials accurately. Our method utilizes sequential ℓ = 0 components to capture isotropic parts that are associated with quantities such as the total permittivity. The sequential ℓ = 2 components are employed to capture anisotropic parts associated with quantities such as birefringence and direction-dependent refractive index. We demonstrate the capabilities of TSENN by considering the frequency-dependent dielectric tensors of a dataset of 1,432 nonmagnetic semiconductors calculated via first-principles methods using OpenMX [37–40]. By training TSENN on the aforementioned dataset, we achieve a MAE of 0.127 in the predicted imaginary part of the dielectric tensor. The high prediction accuracy enables reconstruction of the real part of the dielectric tensor via the Kramers-Kronig (K-K) relations, allowing us to compute additional optical properties derived from the spectra. Accurately predicting the full tensor while preserving its symmetry opens new opportunities for designing materials with tailored anisotropic responses, with applications in photovoltaics, optical modulators, and sensing technologies.

Results
2.1 Dataset
With our focus on small bandgap semiconductors in mind, we selected materials from the Materials project [41] with band gaps ranging from 0.3−3 eV (photon wavelengths of 4,100 nm (infrared) to 413 nm (violet)) that are appropriate for optoelectronic applications [42, 43]. To maintain structural stability, we imposed an energy-above-hull threshold of 0.02 eV/atom. For computational efficiency, we limited our selection to materials with fewer than three elements and fewer than ten atoms per unit cell. Additionally, we excluded magnetic materials as well as materials containing elements with f electrons due to the lack of availability of accurate predefined pseudo-potential sets. These constraints yielded a dataset of 1,432 materials, whose composition is summarized in Fig. 1. Importantly, we did not impose symmetry constraints, such as filtering for specific space groups. As a result, our dataset includes anisotropic materials, whose tensor properties are not intrinsically diagonal. For each entry, we first symmetrized the structure by following the procedure of Ref. [44] and performed first-principles calculations using OpenMX. The dielectric tensors were computed using the Kubo formula for photon energies ranging from 0−30 eV, with an increment of 0.01 eV, resulting in 1,432 frequency-dependent dielectric tensors satisfying the symmetry constraints of the crystal systems considered [45–47]. In this work, the complex dielectric tensor refers to the frequency-dependent relative dielectric tensor, expressed as εαβ(ω) = εαβ(ω)+iεαβ(ω), where εαβ and εαβ denote the real and imaginary parts, respectively. The superscripts αβ index the directional dependence of the spectra. For simplicity, we omit the explicit frequency dependence (ω) hereafter.

To predict the dielectric tensor using an equivariant neural network, we express it in the spherical–harmonics basis, which naturally encodes rotational symmetry. The dielectric tensor is symmetric under the permutation of Cartesian indices, εij = εji, and therefore admits the decomposition
εij = Tαβ Yα Yβ = ε(0) ⊕ ε(2) = Σℓ∈{0,2} Σm=−ℓ^ℓ εm(ℓ) Ym(ℓ), (1)
which separates the scalar trace component (ℓ = 0) from the traceless symmetric tensor component (ℓ = 2). Details of the decomposition procedure are provided in SI Sec. 6. This approach transforms the tensors from the Cartesian basis into the spherical–harmonics basis, making it easier to work with tensor properties in systems with rotational and inversion symmetries. In Eq. (1), note that a rank-2 Cartesian tensor can be decomposed into ℓ = 0 and ℓ = 2 channels, which involve 2ℓ+1 different magnetic quantum numbers m ∈ {−ℓ,...,ℓ} that form a spherical tensor of order ℓ. Since we are dealing with frequency-dependent tensors, we take the Cartesian tensor at each photon energy, perform this change-of-basis operation, and obtain the frequency-dependent spherical harmonic coefficients as our target as shown in Fig. 2. The change of basis and the evaluation of the Wigner 3j coefficients were carried out using the e3nn package [21].

Fig. 2 illustrates the spherical–harmonics decomposition of the dielectric tensor for two representative materials from our dataset, which are shown both in the Cartesian and spherical–harmonics basis. It can be shown through a crystallographic analysis that there are five distinct cases of non-vanishing dielectric tensor components across the seven crystal systems. We highlight two examples. (i) Triclinic compound Br2PdSe6, which lacks any symmetry operations, so that all components of the dielectric tensor are non-vanishing. Consequently, the decomposition of its Cartesian dielectric tensor yields non-zero contributions across all irreducible spherical components, including the scalar channel Y0 and the five rank-2 spherical harmonic components Y2 for m∈{−2,−1,0,1,2}. And (ii) the hexagonal compound AsBaLi in which the higher crystal symmetry imposes constraints on the tensor. Specifically, the sixfold rotational symmetry enforces the presence of diagonal components only with the relationships εxx = εyy ≠ εzz. The xx and yy components (blue and red curves in Fig. 2) are visually indistinguishable within numerical accuracy. All off-diagonal components vanish, consistent with the symmetry constraints summarized in the inset of Fig. 2. As a result, the spherical–harmonics decomposition becomes sparse, with only the Y0 and Y2 components remaining nonzero. Origin of these spherical channels can be traced back to the specific linear combinations of Cartesian components: the scalar component Y0 arises from the trace xx+yy+zz, while the rank-2 spherical harmonic components correspond to the combination −1(xx+yy)+zz [48].

2.2 Equivariant Graph Neural Network Architecture
Our TSENN model extends the e3nn and GNNOpt frameworks, both of which have demonstrated strong performance in sequential prediction tasks such as electron and phonon densities of states, as well as the traces of optical spectra [29, 30, 49]. While TSENN is built on the foundational principles of equivariant GNNs as implemented in e3nn, it incorporates key extensions that distinguish it from GNNOpt. In particular, TSENN leverages higher-order spherical harmonic components to better capture the anisotropy of tensorial properties (Fig. 3). The effectiveness of the output representation, model architecture, and comparison with state-of-the-art methods is further detailed in the Supplementary Information (SI). Accordingly, we construct a GNN that represents each material as a crystal graph, where the atomic and structural information is encoded through node and edge features. Each atomic node i is initially represented by a one-hot encoding of its atomic species z = onehot(atomic number) ∈ R^118, corresponding to the 118 known chemical elements [50]. These atomic types are embedded via a multilayer perceptron to produce initial node-level features fi.

TSENN enforces equivariance to 3D Euclidean transformations by adopting the formalism of tensor field networks [22]. Translational equivariance is achieved through the use of relative atomic positions in message passing, while rotational equivariance is enforced by constraining intermediate features to lie within irreducible representations of the rotation group SO(3). These are indexed by rotation order ℓ ∈ {0,1,2,...} and parity p ∈ {−1,1}. Each feature tensor takes the form V(ℓ,p), where c indexes feature channels and m ∈ [−ℓ,ℓ] indexes spherical harmonic components. Thus, a representation of degree ℓ occupies a tensor space of dimension RNchannels×(2ℓ+1).

These representations are combined using equivariant tensor product operations ⊗ to construct the pointwise convolution layers. This operation can be applied iteratively L times, following the Clebsch–Gordan decomposition as implemented in e3nn:
(U(ℓ1,p1) ⊗ V(ℓ2,p2))(ℓ0,p0) = Σℓ1 Σℓ2 C(ℓ0,m0) U(ℓ1,p1) V(ℓ2,p2), (2)
where ℓ0 ∈ [|ℓ1−ℓ2|, ℓ1+ℓ2] and p0 = p1p2. For brevity, we omit the indices ℓ, p, c, and m when referring to these tensors hereafter.

After the final pointwise convolution layer, we apply mean pooling over all nodes to obtain a global representation of the material. The network output is then partitioned into two parts: a scalar component corresponding to ℓ = 0 and a tensor component corresponding to ℓ = 2. To train the model, we define a composite loss function L, which separately evaluates the scalar and tensor components of the predicted optical spectra. Let N denote the number of materials and n the number of frequency points. For each material entry i, yˆ(i)(ω) and yˆ(i)(ω,m) denote the predicted ℓ = 0 and ℓ = 2 components, respectively, and y(i)(ω) and y(i)(ω,m) denote the corresponding ground-truth values. Here, the hat symbol denotes the predicted tensor, unless explicitly specified otherwise, distinguishing it from the ground-truth tensor obtained via the Kubo formula. Since there are two objectives (L (ℓ = 0) and L (ℓ = 2)) with different magnitudes, we adopt the loss balancing approach for multi-objective training, which allows the network to determine the optimal weighting between objectives dynamically [51]. The effectiveness of this loss balancing is discussed in SI Sec. 2.2. The total loss is defined as:
L = 1/(2σ1^2) L(ℓ=0) + 1/(2σ2^2) L(ℓ=2) + logσ1 + logσ2 (3)
where σ1 and σ2 are learnable parameters that represent the relative uncertainty of each task, with component-wise definitions:
L(ℓ=0) = 1/(Nn) Σi=1^N Σω |yˆ(i)(ω) − y(i)(ω)|, (4)
L(ℓ=2) = 1/(Nn) Σi=1^N Σω Σm=−2^2 |yˆ(i)(ω,m) − y(i)(ω,m)|, (5)
where Mℓ=2 = 2·2+1 = 5 is the number of components for the rank-2 irreducible representation. The model used for prediction was selected based on the best performance according to the validation loss. After obtaining the predicted spherical tensor components from the best model, we transformed them back into the Cartesian tensor representation.

To evaluate the model’s performance, the most direct measure is the mean absolute error (MAE) of the dielectric tensor,
MAE = 1/(6N) Σω Σα≤β |εˆαβ(ω) − εαβ(ω)|. (6)
This ensures that peak position, width, and sign are faithfully captured in the evaluation, without cancellation across frequencies or components. In addition, we compute the per–component errors MAEαβ, as well as a Normalized MAE (NMAEαβ), defined by normalizing MAEαβ to the characteristic magnitude of each tensor component. This scale–aware metric highlights how the prediction errors compare to the intrinsic strength of the typical dielectric response, thereby providing additional physical interpretability.

To isolate anisotropic behavior, we remove the isotropic contribution from the tensor,
εaniso(ω) = ε(ω) − 1/3 Tr[ε(ω)]I. (7)
This decomposition yields the component–wise anisotropic contributions εαβaniso(ω) and their overall magnitudes ∥εaniso(ω)∥F. From these, we define per–component errors MAEαβaniso, the related component–averaged value MAEaniso, and the Frobenius–norm errors MAEnormaniso, to quantify, errors in the overall anisotropy strength. For completeness, we also report additional shape-sensitive metrics such as the first–derivative MAE of individual components MAE′,αβ and the KL divergence for the diagonal components. Detailed definitions of these metrics are provided in the SI Sec. 1.

2.3 Prediction Results
Table 1 summarizes error metrics for dielectric tensor predictions across chemically and structurally diverse test systems. The overall tensor MAE is 0.127 (0.097 median), establishing TSENN’s robust baseline accuracy. When decomposed into isotropic and anisotropic contributions, the averaged anisotropic error MAEaniso is 0.096 (0.074 median), demonstrating reliable capture of directional variations. By contrast, the Frobenius–norm anisotropy error MAEnormaniso is larger at 0.141 (0.111 median), reflecting the stricter requirement of reproducing the total anisotropy magnitude rather than individual Cartesian components. We emphasize that this norm-based measure is the most physically meaningful measure, since it is rotationally-invariant and directly corresponds to the ℓ=2 spherical–harmonics channel. This invariance complements the model’s equivariance, which guarantees the correct transformation of anisotropic components as a rank-2 tensor under rotations.

At the component level, the diagonal elements (xx, yy, zz) exhibit higher errors (0.237–0.259, with NMAE ∼3%) compared to the off-diagonal terms (xy, xz, yz) (0.082–0.105, with NMAE ∼5–7%). This hierarchy is physically natural: the diagonal responses are much stronger in magnitude and receive contributions from both the isotropic (ℓ = 0) and anisotropic (ℓ = 2) channels, whereas the off-diagonal terms arise purely from the ℓ=2 sector and vanish in many high-symmetry systems. Note that we do not apply per-component normalization or weighting: equivariance permits only blockwise scalings that are constant on each SO(3) irreducible subspace. Any differential weighting inside the ℓ=2 block (e.g., different weights for m=±2,±1,0 or per–Cartesian rescaling) would break rotational equivariance.

Beyond absolute and relative magnitudes, values of first-derivative MAEs (0.094–0.105 for diagonal and 0.034–0.039 for off-diagonal terms) and KL divergences (0.106–0.161 for diagonal term) confirm that the model reproduces not only overall intensities but also the line shapes and peak positions. The close agreement between the mean and median values across all metrics further indicates that performance is uniform across diverse systems, without being skewed by outliers. These results demonstrate that TSENN achieves high fidelity in predicting the full dielectric tensor, capturing both the dominant diagonal responses and the more subtle anisotropic features with accuracy and physical consistency.

To visualize the predictions, we stratified the test set based on MAE values and selected representative materials from each quartile, excluding cubic lattices since their isotropic dielectric response is already well understood. This sampling highlights systems with nontrivial anisotropies [Fig. 4(a)]. With the accurate predicted imaginary parts, we used K-K relations to obtain the real parts of the dielectric tensors. Since K-K relations are model-independent, they provide a rigorous test of the internal consistency of the predicted spectra. In Fig. 4(b), we compare the real parts obtained via K-K relations from the predicted imaginary spectra against the reference real parts computed directly using the Kubo formula. An exemplar monoclinic system in the test set is used (Se2Te4), where due to symmetry constraints, only four components of the dielectric tensor εˆij are non-zero: εˆxx, εˆyy, εˆzz, and εˆxz, while the off-diagonal components εˆxy and εˆyz vanish. As shown in the figure, the reconstructed components exhibit excellent agreement with the Kubo-derived results, demonstrating that both the off-diagonal elements (e.g., εˆxz) and the anisotropic diagonal elements (εˆxx, εˆyy, εˆzz) are faithfully reproduced. The other off-diagonal components, which are not shown, are found to be negligible in keeping with symmetry constraints. This analysis validates not only the predictive accuracy of the imaginary dielectric spectrum but also the physical plausibility and self-consistency of the full complex tensor. These results demonstrate the model’s potential for downstream applications, including photonic materials discovery, symmetry-aware property screening, and the inverse design of materials with tailored optical anisotropy.

To visualize the model’s performance and its ability to preserve symmetry across different crystal systems, we randomly selected five materials from the test set, each representing a distinct crystal class. These materials are arranged according to increasing symmetry, from triclinic (minimum symmetry) to cubic (maximum symmetry), see Fig. 5. In the Cartesian tensor representation, the predicted dielectric tensor exhibit symmetry-consistent behavior, correctly reproducing the expected nonzero tensor components in agreement with the ground truth.

To further confirm symmetry preservation, we applied a symmetry-based masking procedure to the predicted tensor at each photon energy (see Methods). While equivariance enforces rotational symmetry at the representation level, the network is not explicitly given point-group information. In practice, small numerical deviations or approximate symmetries in the input structures can lead to residual off-diagonal terms. The masking procedure ensures that the predicted spectra conform to the exact crystal symmetry across the entire spectral range. These distinct symmetries are also reflected in the spherical–harmonics decomposition. Each crystal system preserves different spherical channels that depend on its symmetry. This feature is observable in the spherical tensor plots of Fig. 5, where the colormap indicates the magnitude of each spherical–harmonic coefficients, providing a compact, symmetry-adapted view of the underlying Cartesian dielectric tensor. The prominent red region in the Y0 channel corresponds to a strong isotropic contribution, originating from the trace of the tensor, i.e., εˆxx + εˆyy + εˆzz. In contrast, strong features in the Y2 channel (red or blue in Fig. 5) indicate significant anisotropic contributions from the diagonal components (e.g. −1(εˆxx+εˆyy)+εˆzz). A detailed breakdown of these contributions can be found in Ref. [48]. For example, in a hexagonal material, such as Br6Cd2Cs2, a visible peak appears in the Y2 channel, highlighted in red and blue in Fig. 5, which reflects the anisotropic nature of the diagonal components due to the symmetry condition (εˆxx = εˆyy ≠ εˆzz). This illustrates how the spherical–harmonics basis captures not only the symmetry but also the complexity of the directional responses.

As already pointed out, in the spherical tensor representation, the dielectric tensor can be decomposed into ℓ = 0 and ℓ = 2 channels. Since there are no additional irreducible components at this order, the ℓ = 0 channel remains decoupled as it cannot mix with higher-order terms. Improper mixing between the channels could introduce unphysical features in the predicted spectra. However, because both the input crystal structure and the corresponding dielectric tensor are inherently symmetry-consistent, the model preserves the underlying symmetry constraints during training. As a result, it outputs only the physically allowed tensor components, maintaining the correct form of the spherical decomposition. This symmetry-preserving behavior ensures that no spurious mixing occurs between channels and that the reconstructed Cartesian tensor contains only the symmetry-permitted elements. Consequently, the model captures the intrinsic symmetry of materials without violating crystallographic constraints.

Methods
4.1 Ab initio Calculations
DFT calculations were performed using OpenMX [37–40], which employs numerical pseudo-atomic orbitals as the basis set to expand the one-particle Kohn-Sham wave functions. The high-throughput calculations involved three main steps: generating input files, submitting jobs, and monitoring their progress. The input files were generated using the ASE package, which provides a predefined lookup table for optimized basis selection. However, this table has not been updated since 2013, leading to missing entries for some elements with heavy f-orbital character. We used the PBE-GGA functional in all calculations. We use a custom workflow for job submission and monitoring, implemented with Atomate [55]. Self-consistent field (SCF) convergence was ensured to an accuracy of 3.67×10−8 Hartree or 10−6 eV. All other settings used the default OpenMX parameters. The SCF k-point grid was generated using the pymatgen automatic density function with a target density of 2000 k-points per atom (kppa). This method selects Γ-centered meshes for hexagonal and face-centered lattices, and Monkhorst–Pack grids for all other crystal systems, ensuring consistent and resolution-adaptive Brillouin zone sampling.

After achieving SCF convergence with OpenMX, we extracted the Hamiltonian, overlap matrices, and position operators to evaluate the tensorial relative permittivity. The dielectric tensor was computed within the IPA using the Kubo formula:
εαβ(ω) = εαβ(ω) = χαβ(ω) + δαβ, χαβ(ω) = e^2/ℏ ∫ d^3k/(2π)^3 Σn,m fmn r^nαm r^mβn / (ωmn − ω − iη) (8)
In this work, we denote the real and imaginary parts of the εαβ as εαβ1 and εαβ2, respectively. Here ℏωmn = En − Em denotes the interband energy difference, fmn = f(En) − f(Em) is the difference of Fermi–Dirac occupations, and rαnm is the Berry connection, defined as rαnm = i⟨un|∂α|um⟩, with ∂α = ∂/∂kα the derivative with respect to the Cartesian momentum coordinate α ∈ {x,y,z}, where the Hamiltonian and position operator are projected onto a localized basis and the Brillouin-zone integrals are evaluated efficiently [56–60]. A uniform k-point mesh was constructed with reciprocal-space resolution of 0.05 Å−1 along each reciprocal lattice vector, ensuring consistent Brillouin-zone sampling across all materials. The smearing parameter η, introduced to improve numerical convergence of the optical response, was fixed at 0.05 eV in all calculations. This value corresponds to typical carrier relaxation times (∼10 femtoseconds) in semiconductors [61–63]. This procedure is equivalent to constructing a Wannier-based tight-binding representation of the Hamiltonian, which provides a localized framework for efficient evaluation of the dielectric response. The present calculations are performed within the IPA, i.e. based on Kohn–Sham eigenvalues and wave functions without including quasiparticle (GW) corrections or electron–hole interaction (excitonic) effects [64]. Consequently, discrepancies between the computed dielectric spectra and the corresponding experimental results are expected.

4.2 Machine Learning
The spectra were computed with a resolution of 0.01 eV, with a total of 3,001 points. Including all points in the training process would significantly increase computational cost. To make training more efficient, we performed downsampling by interpolating the data along the photon energy grid to 300 points, which still captures the spectral response reasonably. We used an 80:10:10 split for training, validation, and testing datasets. Given our focus on crystal symmetries, we implemented balanced loading by analyzing the distribution of crystal systems across all splits to ensure even representation, see SI Sec. 7 for details.

Group theoretical analysis reveals five distinct cases based on crystalline symmetries. This enables us to construct a mask M that identifies the symmetry-allowed, non-vanishing components of the tensorial spectra. The mask is derived by initializing a 3×3 matrix of ones, 1_3×3, and applying the set of rotation operations O = {R1, R2, ..., RNR} associated with the space group symmetry of the crystal. The final mask is computed by averaging over the transformed matrices:
M = 1/NR ΣRi∈O Ri⊤ 1_3×3 Ri (9)
For example, in a monoclinic system, this process yields the following mask:
(1 0 1)
M = (0 1 0) (10)
(1 0 1)
This mask highlights the non-vanishing tensor components, such as xx, yy, zz, and xz, as dictated by the symmetry constraints. To identify any symmetry violations, we define the residual mask as
Residual mask = ε(1_3×3 − M) = (0 ε 0)
(ε 0 ε) (11)
(0 ε 0)
where ε = 10−6 serves as the detection threshold. We apply this residual mask across the photon energy grid to verify the input tensors before training and the predicted tensors after training. Our results confirm that all tensors respect the imposed symmetry constraints, demonstrating that the model maintains equivariance throughout training without symmetry breaking.

We performed a comprehensive hyperparameter search using Optuna [65] to optimize performance of the model. The multiplicity of irreducible representations m used in the convolution layers was varied from 16 to 64, with 32 selected as the optimal value. The number of pointwise convolution layers L was tuned between 1 and 4, with the final model using 4. The embedding feature vector length was explored from 16 to 160, with 128 yielding the best performance. Dropout probability was tested at 0, 0.1, 0.2, and 0.4 with 0.4 selected as the optimal choice. Batch sizes from 1 to 16 were systematically evaluated, and a size of 8 was selected as the best tradeoff between computational efficiency and training stability. Smaller batches (e.g., size 1) led to unstable optimization due to insufficient statistical averaging per update. For very large batches (>32), the anisotropy was found to become less pronounced compared to smaller ones, likely because the optimization tends to get trapped in local minima and averages out sharper anisotropic signals through stochastic smoothing. The hyperparameters of the AdamW optimizer were carefully tuned. Learning rates between 5×10−4 and 5×10−2 were explored, with 1×10−2 achieving the lowest validation loss. To further enhance convergence, we applied PyTorch’s CosineAnnealingWarmRestarts scheduler, which gradually reduces the learning rate according to a cosine decay schedule. The initial number of iterations before the first restart was set to T0 = 10, with subsequent settings kept at their default values. This smooth, non-monotonic decay enables larger updates during the early training stages and finer adjustments at later stages, helping the model escape shallow minima and converge more robustly [66].

To contrast the computational efficiency, we benchmarked our workflow against a representative ab initio calculation. For this purpose, we calculated the dielectric tensor of Ag3AsS4 (Materials Project ID: mp-9538, 16 atoms per unit cell). The DFT self-consistent calculation converged in 119 s using a hybrid MPI/OpenMP setup on 2 AMD EPYC 7763 CPUs (256 threads in total) with default workflow settings. The subsequent dielectric tensor evaluation, carried out with Brillouin-zone integration parallelized over 128 processes, required 22 min 9 s. Including workflow orchestration overhead, the full pipeline–from job launch to data storage–took 26 min 32 s. By contrast, the GNN-based workflow was orders of magnitude faster: graph construction required 6.96±0.01 ms, forward inference 695.20±0.01 ms, and K-K relations took 570.04±0.01 ms, giving a total runtime of 1.18 s on a single NVIDIA RTX 4090 GPU. In other words, our approach reduces a half-hour DFT calculation to just over one second.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript addresses the prediction of all tensor components of spectra, i.e., frequency-dependent quantities such as the dielectric function, using equivariant GNNs (Graph Neural Network).\nSeveral groups have recently shown that GNNs can predict spectra in principle.\nIn parallel, several groups have shown that tensorial properties can also be predicted well.\nThis leaves a certain gap, namely the prediction of frequency-dependent tensorial properties.\nThe present work more or less fills this gap.\nThe relevant community will not be surprised by this closure of the gap, but will welcome it, so that publication of the results is in principle to be welcomed.\nThe qualification “in principle” has several aspects.\nPerhaps the most important is the question of whether Nature Communications is the right journal, since relatively specialized prior knowledge is required to recognize the described gap as such and to appreciate its closure.\nFurthermore, I also find, not least because of the choice of presentation of the results, that the successful description of the tensor character is not as convincing as one would wish.\nIn addition, there are a surprising number of careless errors or presentation weaknesses, which are not decisive in terms of content but are annoying.\n---\n*Major issues*\n(1) As stated above, GNN have proven to be able to predict scalar spectra such as the traces of permittivity tensors.\nTherefore, the central criterion for evaluating the present manuscript must be whether the non-diagonal components such as chi..xz [.. marks subscripts henceforth] are also well predicted.\nStrictly speaking, there is only one place in the whole manuscript where useful quantitative information can be found,\nnamely in the first two rows of the two left columns in Fig. 5.\nThe specific question is whether, when viewed impartially, one would say that the green lines of the corresponding panels match well or not.\nUnfortunately, I would answer this question in the negative!\n(2) This also highlights another problem with the representations throughout the manuscript: the off-diagonal components,\nwhich are being investigated for the first time, are almost invisible in view of the much larger diagonal components.\n(3) The authors provide a measure of error/uncertainty, namely the quantity overline(overline(chi)) in Eq. (10).\nThis analysis is flawed for at least ,,, reasons.\n(i) The chi-values are summed of frequencies:\nRoughly speaking, the dominant feature of the permittivity of all considered materials is a more or less isotropic peak.\nIf the position and width and shape of this peak were predicted completely *wrong*, the defined error measures\nMAE = |overline(overline(chi))^prediction - overline(overline(chi))^target|\ni.e. (Eq. (11)) and the corresponding 'relative error' (Eq. (12)) would still be zero, as long as the overall weight (~ oscillator strength) is predicted correctly.\n(ii) Positive and negative off-diagonal elements are summed:\nIn contrast to the diagonal elements, off-diagonal elements can have either sign.\nE.g. chi..xz will change sign, if we replace z by -z.\nThis purely mathematical transformation changes the value of the quantity defined in Eq. (11) and thus the error estimate.\n(iii) For many uniaxial materials, the leading deviation from isotropy defined as trace(chi)/3 are finite values of the quantity:\nchi..zz - 0.5*(chi..xx+chi..yy).\nThe considered quantity overline(overline(chi)) is complete blind with respect to this anisotropy.\nSummary of (i), (ii) and (iii):\nOne should consider\nsum-over-omega sum over components of absolute value of ( prediction - target )\ninstead of\nabsolute value of sum-over-omega sum over components of ( prediction - target )\n(iv) Most considered components do not have off-diagonal tensor elements:\nTherefore, these error measures are dominated by compounds where we know that GNN can be successful.\nSince the new contribution of the manuscript are the anisotropic components, it is not useful to present statistics which are\ndominated by compounds not showing such components.\n(v) Even for compounds that show anisotropies as quantified by e.g. chi..zz-0.5*(chi..xx+chi..yy) and chi..xz, these are\ndominated in the defined error measures by the much larger isotropic contribution trace(chi)/3.\nSummary of (iv) and (v):\nOne should give results for, e.g. the quantities chi..zz-0.5*(chi..xx+chi..yy) and chi..xz including only those compounds of Fig.\n1a, where symmetry does not force those quantities to vanish.\nI am quite convinced that if someone had accidentally drawn Figure 4 using the isotropic mean value and completely\ndisregarding all off-diagonal components, no one would have noticed.\nThis implies that Figure 4 has no bearing on the subject of the paper.\n(4) I get the impression that the authors consider it a success that anisotropic and in particular off-diagonal components\nvanish for the more symmetric compounds.\nIn this context, I have two questions going in different directions:\n(4a) Isn't it obvious that predictions of an equivariant network preserve the symmetry.\nIn other words: Isn't it trivial that anisotropic and in particular off-diagonal components vanish for the more symmetric\ncompounds? Or am I mistaken?\n(4b) I like the comparison behind Supplementary Figure 1.\n(5) Some early papers showing the ability of GNN to predict scalar spectra have not been cited, e.g.\nM. Grunert et al, Phys. Rev. Materials 8 (2024), L122201; A. Ibrahim and C. Ataca, ACS Appl. Materials & Interfaces 16\n(2024), 41145-41156\n(6) Why is section 4.2 called 'tight binding calculations'?\nI see nothing in the entire section that has anything to do with tight binding except for the general references [41,49-52] and a\ncertain similarity between the Wannier orbitals mentioned there and the tight-binding picture.\n(7) Many people, including https://en.wikipedia.org/wiki/Permittivity\ninterpret 'permittivity' as epsilon(omega) and write\nepsilon = epsilon..0 * (1+chi).\nThe present manuscript deviates from this usage in two ways:\n(i) chi is called permittivity\n(ii) chi includes epsilon..0 (see Eq. (13)).\nThese choices should be stated explicitly in order to avoid confusion.\nEven better would be to present the results for dimensionless quantities, such es epsilon..rel = epsilon / epsilon..0 or the\nconventionally defined chi.\nMost people have a feeling for dimensionless epsilon..rel values, but not for values in units of milifarad/m.\n(8) The manuscript discusses only spectra calculated within the Kohn-Sham independent particle picture.\nThe fact that and why these might differ dramatically from experimental spectra should at least be mentioned.\nIn SUMMARY,\nI can not recommend publication of these results in Nature Communications,\nbut I encourage the authors to publish their results with much more focus on the an-isotropic and off-diagonal tensor\ncomponents.\nIn particular, I strongly recommend to use error measures that come closer to what humans would do if the look at Fig. 5 to\ndecide whether they consider the green lines to be similiar or not.\n---\n*Exemplary, not exhaustive list of minor issues*\nThe subsection \"Author contribution\" states \"T.-W. H., Z.F. and Q.Y. conceived the experiments.\"\nGiven that no experiments are presented, this surely is a cut\\&paste error.\nThe reference list is full of names of persons or elements which are not capitalized: \"weyl metal\", \"van der waals\", \"from h to\nkr\".\nThis can easily be fixed, but should have been fixed before submission.\nSimilarly: If the 'S' in Ref. [54] \"Sgdr: Stochastic.. \" is capitalized, one should capitalize the Name of the Method 'SGDR' as\nwell.\nIn the abstract, the authors state \"our model achieves a mean absolute error (MAE) of 21.181 millifarads per meter (mF/m)\".\nWe teach our students that specifying an error estimate with 10^(-5) accuracy makes no sense at all.\nSince I have just quoted this sentence, I would like state again that it would probably make more sense to consider\ndimensionless quantities. Honestly, how many readers are familiar enough with permittivity values to appreciate this\nnumerical error value? In contrast, many know that epsilon..relativ values are of order 1-10.\nFurthermore, in an abstract in a renowned journal such as Nature Comm., it is not necessary to explain that ‘mF’ stands for\n‘millifarad’ in the context of permittivity.\nIn addition to specific references to concrete points, this work naturally also contains more general references.\nHowever, some of these seem to have been chosen extremely arbitrarily:\nFor example, the introduction begins with \"... which together form the foundation for interpreting a wide range of optical\nphenomena in solid-state systems [1, 2].\"\nThe first reference \"[1] Lucarini, V. Kramers-Kronig relations in optical materials research (2005)\" refers to a book on a very\nspecial aspect of optical spectra, without indicating the publisher.\nThe second reference \"[2] Boyd, ...Nonlinear Optics, 1097–1110 (Springer..., 2023)\" refers to a chapter of a book, specifying\nthe page numbers and the name of the chapter, but not the name of the book.\nAnd why is a chapter on 'nonlinear' optics cited if the whole manuscript addresses 'linear' optics?\nI find the references for \"we selected materials ... appropriate for optoelectronic applications [34, 35]\"\nHow did ref. [34] \"Band parameters for III–V compound semiconductors and their alloys\" influence the choice of the 1432\nbulk semiconductors considered in this study.\nFig. 6(b) includes as inset a linear fit of chi..xz vs. shear-strain angle.\nNaively, I would aspect a sine relation.\nIf I look at the data, I get the feeling that a sine shape starting slightly deeper and then leveling off slightly could possibly fit\neven better.\nIt should also be relatively simple for the authors to compute the actual values of the sheared/stretched crystals, to verify\nwhether the model predicts the correct trend.\nGiven that an equivariant model should by definition only predict the correct tensor elements, the purpose of the section on\nmasking the resulting tensor is not clear to me.\nThe prediction time of TSENN should also be contrasted with the time taken by the equivalent ab initio calculation.\nSome minor questions regarding the choice of different parameters:\nWhy are the chosen batch sizes of the model so small - is the VRAM usage of the model that high?\nIf so, that should be indicated as well\nWhy was the broadening of the dielectric function chosen to be 50 meV, a relatively small value?\n(Remarks on code availability)\nWe looked only superficial at the code and the parameter choices.\nWhile some code has been made openly available, the code for the ablation experiments is missing.\nWhile this is very important for possible users of TSENN, it would be beneficial for the review as well.\nFurthermore, the workflows for the generation of the optical data appears to be missing too.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.\n(Remarks on code availability)", "reviewer_3": "Hsu et al. have demonstrated a powerful neural network architecture, namely TSENN, that aims to relate the crystal structure\nto the full anisotropic frequency-dependent dielectric permittivity. The demonstrated architecture shows high prediction\naccuracy, confirmed by the KK relations with capabilities of maintaining the anisotropy. Furthermore, authors have presented\ntwo case studies of strain-induced symmetry breaking and the resulting change in dielectric permittivity.\nThe data-driven results are timely and interesting, and particularly useful for materials screening for both new materials\ndiscovery as well as application purposes – a helpful methodology for both academia and industry. However, the advantage\nof simplifying the symmetry processing and accuracy in permittivities, particularly below/near bandgap, would require further\nclarification. Detailed questions are shown below:\n1. What are the advantages of using spherical harmonic decomposition compared with using dielectric permittivity? For non-\nmagnetic materials, both approaches yield the same number of variables for all seven lattices, given symmetric off-diagonal\ndielectric components (i.e., x_ij = x_ji) and rotational symmetry. Thus, to me, spherical harmonic decomposition does not\noffer significant differences in the anisotropy results and improvement in reducing the computational complexity. For\nexample, a uniaxial system has 2 unique dielectric permittivities (xx, zz), 4 for monoclinic (xx, xz, yy, zz), and 6 for triclinic,\nand the same variables are generated from spherical harmonic decomposition. If no significant advantages, what’s new\nregarding this approach?\n2. Typically, the dielectric permittivity below the plasma frequency matters more as that highlights distinct material\ncharacteristics. Looking at figure 4b, the spectrum shows more deviation in the absorption spectrum, while better agreement\nabove plasma frequency (where fewer material distinctions will be shown).\na. Taking more frequencies above the plasma frequency reduces the statistical weight of the absorption spectrum; thus, I\nthink authors should have a statistical discussion and figure focusing on frequencies below the plasma frequency.\nb. Accurately predicting the bandgap has been a long-standing topic in the community. It would be beneficial if authors could\ncompare the bandgap from predicted dielectric permittivities (e.g., Tauc method[1]) vs the calculated bandgap to highlight\nthe prediction accuracy.\n3. Magnetic materials can have non-symmetrical dielectric permittivity. Would current methods capture all 9 dielectric\npermitivities without assuming symmetrical components?\n[1] B. D. Viezbicke, S. Patel, B. E. Davis, D. P. Birnie III, Phys. Status Solidi B 2015, 252, 1700.\n(Remarks on code availability)", "reviewer_4": "This paper presents a model that predicts the full, frequency-dependent dielectric tensor directly from structure using an\nequivariant GNN (E(3)NN). The outputs are decomposed into ℓ=0 (isotropic) and ℓ=2 (anisotropic) spherical components.\nThe dataset includes 1,432 non-magnetic semiconductors from the Materials Project, recalculated with OpenMX (numerical\npseudo-atomic orbitals). The imaginary part of the dielectric tensor is computed via the Kubo–Greenwood formalism for 0–30\neV (0.01 eV step, but downsampled by interpolating the data onto a 300-point photon energy grid for GNN training). The\nauthors report a test MAE of 21.181 mF/m and use Kramers–Kronig to get the real part. They also show symmetry checks\nunder strain/shear that indicate the model respects crystallographic constraints. I suggest publication of this manuscript after\nminor revision. My comments are given below:\n1. In the intro, the authors write:\n“Notably, existing studies mostly focus exclusively on the dielectric constant, i.e., the ω → 0 static limit.\nIn contrast, the joint prediction of tensorial structure and spectral behavior—specifically, the full frequency-dependent dielectric tensor—remains a largely unexplored area. In this work, we address this gap by leveraging the expressive power of higher-order components within the equivariant GNNs to directly predict the full tensorial, frequency-dependent permittivity, enabling a more comprehensive understanding of the anisotropic optical response of crystalline materials.”\nI think this part would benefit from a quick mention of the literature. It might be fair to say that for tensorial spectra, the area\nhas been explored before but more indirectly and in a fundamental way, for example in equivariant machine-learning\nHamiltonians (https://www.nature.com/articles/s41524-023-01130-4, https://www.nature.com/articles/s41467-023-38468-8),\nbut not specifically for the frequency-dependent dielectric spectrum. There *are* already some papers that address crystal\nstructure-to-spectrum predictions for the frequency-dependent dielectric function, absorption spectra, or solar absorbance\nusing equivariant or invariant GNNs (https://pubs.acs.org/doi/full/10.1021/acsami.4c07328,\nhttps://advanced.onlinelibrary.wiley.com/doi/full/10.1002/adma.202409175). I think these deserve a brief mention up front,\nboth to acknowledge the previous work and to clarify that your contribution is the full tensorial extension using an equivariant\nGNN, just to sharpen the novelty claim.\n2. MAE alone as an error metric can miss important spectral details. Two predictions might have the same MAE, but differ a\nlot in terms of peak positions, heights, band gap onset, or static-limit values. Please consider adding one or two spectrum\nshape-aware metrics on the test set (for example, first-derivative MAE, or a shape distance like earth-mover distance or KL\ndivergence on a PDF-normalized spectrum, or other relevant metrics) and report these in a small table or figure. This would\nback up validation about spectral fidelity.\n3. Since the outputs are split into ℓ=0/ℓ=2, it might be helpful to see/analyze an ablation where you train only on ℓ=0 and then\nonly on ℓ=2. This may show how much accuracy comes from the anisotropy vs. isotropy channels, and whether there’s any\nleakage between channels.\n(Remarks on code availability)\nMy comments are reconsiled and given in \"comments for authors\"", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.\n(Remarks on code availability)\nThe TSENN code is accompanied by a clear README, includes necessary scripts and notebooks, and appears to be a\nwell-organized, reproducible resource for the community, although I have not personally run it."}, "author_response": "Letter to the Reviewers\nDear Reviewers,\nWe would like to begin by sincerely thanking you for the time and effort you have invested in reviewing\nand handling our manuscript. We greatly appreciate the three detailed and knowledgeable reports, which\nraised important points ranging from technical clarifications in figures and captions to broader\nsuggestions for strengthening the framing and future directions of this research. Your constructive\nfeedback has been invaluable in helping us substantially improve the clarity, rigor, and overall\npresentation of the manuscript.\nIn the following letters, we address all comments in point-by-point detail: reviewer comments are\nreproduced in black text, our detailed responses are provided in red text, and the corresponding\nrevisions are outlined in blue, with line number referring to the compiled PDF and highlighted in red\nwithin the manuscript for easy reference. In addition to implementing the requested changes, we also\ntook the opportunity to correct minor typographical errors and to standardize references to further\nimprove the readability and consistency of the manuscript.\nReviewer Comments:\nReviewer #1 (Remarks to the Author):\nThe manuscript addresses the prediction of all tensor components of spectra, i.e., frequency-dependent\nquantities such as the dielectric function, using equivariant GNNs (Graph Neural Network).\nSeveral groups have recently shown that GNNs can predict spectra in principle.\nIn parallel, several groups have shown that tensorial properties can also be predicted well.\nThis leaves a certain gap, namely the prediction of frequency-dependent tensorial properties.\nThe present work more or less fills this gap.\nThe relevant community will not be surprised by this closure of the gap, but will welcome it, so that\npublication of the results is in principle to be welcomed.\nThe qualification “in principle” has several aspects.\nPerhaps the most important is the question of whether Nature Communications is the right journal, since\nrelatively specialized prior knowledge is required to recognize the described gap as such and to\nappreciate its closure.\nFurthermore, I also find, not least because of the choice of presentation of the results, that the successful\ndescription of the tensor character is not as convincing as one would wish.\nIn addition, there are a surprising number of careless errors or presentation weaknesses, which are not\ndecisive in terms of content but are annoying.\nAuthor response to Reviewer #1:\nWe sincerely thank the reviewer for the thoughtful and constructive assessment of our work. We\nappreciate the recognition that our study addresses the outstanding gap between spectral property\nprediction and tensorial property prediction—namely, the prediction of frequency-dependent tensorial\nquantities. While the community may not be surprised by this natural extension, we agree with the\nreviewer that providing a concrete demonstration of such capability is both timely and valuable. We\nhope that our work will serve as a foundation for further developments in this direction.\nWe also appreciate the reviewer’s comment regarding the breadth of readership at Nature\nCommunications. While our study necessarily involves concepts from both tensorial spectra and\nmachine learning—particularly equivariant graph neural networks—this dual perspective is central to\nour goal of bridging these two communities. To enhance accessibility, we intentionally focused on a\nrelatively simple yet representative tensorial response (i.e. dielectric function) and adopted one of the\nmost transparent equivariant architectures, accompanied by step-by-step documentation, open data, and\nfully reproducible codes. We believe that a general framework capable of unifying tensorial spectra with\nequivariant machine learning carries broad significance, with potential impact across materials physics,\nchemistry, and engineering. In response to referee’s comment, we have revised the Introduction and\nDiscussion to better motivate the problem and to highlight the relevance of our approach beyond the\nimmediate expert community.\nWith respect to our presentation, we are grateful for the reviewer’s candid feedback. We have carefully\nrevised the manuscript to address the identified weaknesses, corrected careless errors, and improved the\nclarity of all text, figures, and references. We have also expanded the presentation of tensorial results by\nadding discussions of frequency-resolved benchmarking, spherical-harmonic error decomposition, and\nadditional strain validation tests to make our description of anisotropy more convincing and transparent.\nWe hope these substantial revisions address the reviewer’s concerns and improve both the clarity and\naccessibility of the manuscript.\nDetailed responses to reviewer questions are given below:\nQuestion 1. As stated above, GNN have proven to be able to predict scalar spectra such as the traces of\npermittivity tensors. Therefore, the central criterion for evaluating the present manuscript must be\nwhether the non-diagonal components such as 𝜒 are also well predicted. Strictly speaking, there is\n!\"\nonly one place in the whole manuscript where useful quantitative information can be found, namely in\nthe first two rows of the two left columns in Fig. 5. The specific question is whether, when viewed\nimpartially, one would say that the green lines of the corresponding panels match well or not.\nUnfortunately, I would answer this question in the negative!\nWe thank the reviewer for this thoughtful comment and agree that the presentation of off-diagonal\ncomponents in Fig. 5 could be clearer. The intention of this figure was to demonstrate that our model\ncan capture the presence and qualitative behavior of the off-diagonal dielectric components (e.g.,𝑥𝑧)\nwithout any preferential sample selection—the materials shown were drawn directly from the test set.\nThese examples illustrate that the model correctly identifies symmetry-allowed off-diagonal responses\npurely from the geometric information encoded in the input structure.\nTo address the reviewer’s concern, we have now included additional quantitative analyses. When the\nsame monoclinic material (N₄Se₆Y₈ in main text Fig. 5) is stratified by its anisotropy-based error metric\n(MAE$'()), it falls within the third quartile (Q3) of the test-set distribution, indicating moderate\n#$%&'\ndifficulty rather than an outlier, see Response Fig. 1. We also provide new examples of other\nmonoclinic materials spanning all quartiles (Response Fig. 2) to demonstrate consistent performance\nacross anisotropy levels. Furthermore, per-component errors for both diagonal and off-diagonal tensor\nelements are summarized in Response Table 1 (embedded as Table 1 in main text), with a mean\nabsolute error (MAE) of 0.091 and normalized MAE (NMAE) of 7.3% for the off-diagonal channels—\nvalues computed exclusively over symmetry-allowed structures in the test set. The details and the\neffectiveness of these metrics are discussed together in our Response to Question #3 below, to\nwhich we refer the reader in the interest of brevity.\nWe note that predicting off-diagonal elements is intrinsically more challenging because, unlike the\ndiagonal terms, each off-diagonal response arises from a single independent physical channel rather than\nmultiple degenerate pathways. Nonetheless, the quantitative results confirm that our model\nsystematically captures these weaker anisotropic responses, validating its ability to predict both diagonal\nand off-diagonal tensorial spectra in a symmetry-consistent manner.\nImportantly, we cannot arbitrarily adjust the relative weighting of individual tensor channels (e.g. assign\ndifferent weights to ℓ = 2,𝑚 = 2 and ℓ = 2,𝑚 = 1 channels) to amplify the off-diagonal learning\nperformance, as such intervention would explicitly break the rotational equivariance constraint that\nunderpins the model’s physical fidelity. Any preferential weighting will distort the tensor transformation\nproperties and introduce unphysical components, thereby compromising the self-consistency of our\nequivariant formulation.\nResponse Figure 1: The monoclinic system N Se Y in main text Fig. 5 illustrates the comparison between the target (left\n4 6 8\npanel) and predicted (right panel) imaginary parts of dielectric spectra. When stratified according to the 𝑀𝐴𝐸\"%&' metric we\n!\"#$% \ndefined, it falls within the third quartile (Q3) of the test-set distribution, indicating that it presents a moderately difficult case\nfor our model but that it is not an outlier.\nResponse Figure 2: Monoclinic materials spanning all quartiles to demonstrate consistent performance across the anisotropy\nlevels. Per-component errors in both the diagonal and off-diagonal tensor elements are summarized in Response Table 1, with\na mean absolute error (MAE) of 0.091 and normalized MAE (NMAE) of 7.3% for the off-diagonal channels; values were\ncomputed exclusively over the symmetry-allowed structures in the test set.\nSummary of revisions made based on preceding discussion:\n• Main text, Table 1.\n• Main text, line 219—243.\nQuestion 2. This also highlights another problem with the representations throughout the manuscript:\nthe off-diagonal components, which are being investigated for the first time, are almost invisible in view\nof the much larger diagonal components.\nWe appreciate the reviewer’s insightful comment. The apparent invisibility of the off-diagonal\ncomponents arises from the intrinsic nature of the dielectric response function itself, rather than from\nany limitation of our representation. The dielectric tensor is typically dominated by its diagonal terms,\nwhich quantify the polarization along the principal crystal axes, while the off-diagonal terms encode the\nweaker cross-couplings between the orthogonal directions. From our first-principles dataset, the\ncharacteristic magnitudes—defined as the frequency-maximum of each tensor component averaged\nacross all materials—show that the off-diagonal components are on average about five times smaller\nthan the diagonal ones, as summarized below:\nThis intrinsic imbalance makes simultaneous visualization of all six tensor components challenging on a\nsingle scale. Nevertheless, our equivariant model captures both the diagonal and off-diagonal\nanisotropies faithfully—including the symmetry-forbidden vanishing of specific terms in high-symmetry\nsystems (e.g., trigonal, orthorhombic). Compared to non-equivariant baselines, which can only capture\nthe isotropic behavior, our framework reproduces these subtle directional dependencies that reflect the\nunderlying crystal symmetries.\nSummary of revisions made based on the preceding discussion:\n• Added per-component spectra for a representative material, see Fig. 4(b) of the main text.\n(embedded as Response Fig. 3 in Question 3)\n• Additional examples provided in SI Fig. 10 (Response Fig. 4 in Question 3).\nQuestion 3. The authors provide a measure of error/uncertainty, namely the quantity\noverline(overline(chi)) in Eq. (10).\nThis analysis is flawed for at least ,,, reasons.\n(i) The chi-values are summed of frequencies:\nRoughly speaking, the dominant feature of the permittivity of all considered materials is a more or less\nisotropic peak.\nIf the position and width and shape of this peak were predicted completely *wrong*, the defined error\nmeasures\nMAE = |overline(overline(chi))^prediction - overline(overline(chi))^target|\ni.e. (Eq. (11)) and the corresponding 'relative error' (Eq. (12)) would still be zero, as long as the overall\nweight (~ oscillator strength) is predicted correctly.\n(ii) Positive and negative off-diagonal elements are summed:\nIn contrast to the diagonal elements, off-diagonal elements can have either sign.\nE.g. chi..xz will change sign, if we replace z by -z.\nThis purely mathematical transformation changes the value of the quantity defined in Eq. (11) and thus\nthe error estimate.\n(iii) For many uniaxial materials, the leading deviation from isotropy defined as trace(chi)/3 are finite\nvalues of the quantity: chi..zz - 0.5*(chi..xx+chi..yy).\nThe considered quantity overline(overline(chi)) is complete blind with respect to this anisotropy.\nSummary of (i), (ii) and (iii):\nOne should consider\nsum-over-omega sum over components of absolute value of ( prediction - target )\ninstead of\nabsolute value of sum-over-omega sum over components of ( prediction - target )\nWe thank the reviewer for their detailed and constructive comments. We fully agree that the previous\nerror definition (Eqs. (11–12)) assigned a disproportionate weight on isotropic contributions and\nobscured the anisotropic/off-diagonal components that are the focus of this work.\nQuestion on MAE: Following the reviewer’s suggestion, we now compute the mean absolute error\n(MAE) between the predicted and target spectra, i.e.\n6𝑁\n* \n* +,- \nThis ensures that peak positions, widths, and signs of the spectra are faithfully captured in the\nevaluation, without cancellation across frequencies or components.\nSummary of revisions made based on the preceding discussion:\n• Main text, line 202—209.\n(iv) Most considered components do not have off-diagonal tensor elements:\nTherefore, these error measures are dominated by compounds where we know that GNN can be\nsuccessful. Since the new contribution of the manuscript are the anisotropic components, it is not useful\nto present statistics which are dominated by compounds not showing such components.\n(v) Even for compounds that show anisotropies as quantified by e.g. chi..zz-0.5*(chi..xx+chi..yy) and\nchi..xz, these are dominated in the defined error measures by the much larger isotropic contribution\ntrace(chi)/3.\nSummary of (iv) and (v):\nOne should give results for, e.g. the quantities chi..zz-0.5*(chi..xx+chi..yy) and chi..xz including only\nthose compounds of Fig. 1a, where symmetry does not force those quantities to vanish.\nI am quite convinced that if someone had accidentally drawn Figure 4 using the isotropic mean value\nand completely disregarding all off-diagonal components, no one would have noticed.\nThis implies that Figure 4 has no bearing on the subject of the paper.\nQuestion on Figure 4: We acknowledge the original metrics used in old Figure 4 were not optimal as\nthey were shrinking the information too much. We have now introduced a new Figure 4 (Response Fig.\n3) in the main text based on the new MAE introduced above as the evaluation metric. In panel (a), we\nshow a cumulative kernel-density-estimator (KDE) plot of MAE values and stratified sampling across\nthe quartiles (three samples per quartile, 12 in total). Cubic lattices are excluded, since they are already\nwell understood and expected to fit well. We instead focus on materials that exhibit anisotropy. Fig. 4(b)\nprovides a zoomed-in view in an alternative representation, highlighting that accurate predictions of the\nimaginary part also translate into excellent agreement for the real part via the Kramers–Kronig relation.\nNotably, Fig. 4(b) demonstrates that both the off-diagonal (e.g., 𝑥𝑧) and the anisotropic diagonal\ncomponents (𝑥𝑥,𝑦𝑦,𝑧𝑧) are faithfully recovered, with MAE values ranging from 0.082 to 0.232. An\nadditional figure that samples a monoclinic system in the test set is also embedded in the revised SI as\nSupplementary Fig. 10, to present the off-diagonal fidelity of the model. We have moved the old\nFigure 4 and the related discussion to the SI under Supplementary Fig. 2, where we also clarify the\nnuances of the metrics we use.\n(a)\n(b)\nResponse Figure 3 (Figure 4 in main text). (a) Distribution of component-wise MAE values across the test set (left) and\nrepresentative spectra obtained by stratified sampling across quartiles (right). From the cumulative kernel-density-estimator\n(KDE) plot, we randomly select three systems in each quartile (𝑄 −𝑄 ), resulting in 12 representative examples that span\n( ) \nthe full error distribution. Cubic lattices are excluded, as their isotropic spectra is already well studied, and we instead focus\non systems that exhibit finite anisotropy. For each sampled material, the target and predicted tensor components are shown\nside by side, with the corresponding MAE reported in each panel. This stratified sampling strategy ensures that both high-\nand low-error cases are represented, providing a balanced and unbiased view of model performance across crystal systems\nand anisotropy strengths. (b) Real parts of the dielectric tensor components for the monoclinic system Se Te from the test\n4 2 \nset. Black curves show the ab initio dielectric spectra computed within the IPA using the Kubo formula, while red curves\ndenote spectra reconstructed from the predicted imaginary part via the K-K relation. The MAE*+ is reported in each panel,\nwith values ranging from 0.082 to 0.232, demonstrating strong agreement across all non-vanishing components\n(𝜀,,,𝜀--,𝜀..,𝜀,.). Notably, both the off-diagonal element (𝜀,.) and the anisotropic diagonal terms are faithfully recovered.\nResponse Figure 4. (Supplementary Figure 10.) Comparison between the real part of the permittivity tensor obtained via\nK-K relations from the model-predicted imaginary spectra and the direct calculation using the Kubo formula. Excellent\nagreement is observed across all non-vanishing components. Components such as 𝜀,-and 𝜀,. remain zero due to symmetry\nand thus do not contribute to the overall MAE.\nSummary of revisions made based on the preceding discussion:\n• Main text, line 244 – 260.\n• Main text, Fig. 4.\n• SI Fig. 10.\n• SI line 65—85.\n• SI Fig. 2.\nAnisotropy-focused metric (Question (iv)): To isolate the anisotropic behavior, we adopt the\nanisotropy norm introduced in our original manuscript, which was previously discussed in the SI, but it\nis now incorporated in the main text, by removing the isotropic contribution from the tensor,\n1\n𝜀 (𝜔) = 𝜀(𝜔)− Tr[𝜀(𝜔)]𝐼.\naniso 3\n+- \nThis yields the component-wise anisotropic contributions 𝜀 (𝜔), as well as the associated overall\naniso \nmagnitude ‖𝜀 (𝜔)‖ . We thus define per–component errors MAEnorm , the component-averaged\naniso . aniso \nvalue MAE , and the Frobenius–norm error MAEnorm , in order to quantify discrepancies in the\naniso aniso \noverall anisotropy strength.\nThis measure removes the isotropic trace contribution and isolates both the diagonal anisotropy (e.g.\n/\n𝑧𝑧− (𝑥𝑥 +𝑦𝑦)) and the off-diagonal components (e.g. 𝑥𝑧).\nanisotropy across the frequency domain. Note, however, that this metric condenses the tensor into a\nsingle scalar quantity and therefore captures only the magnitude of anisotropy. For this reason, we\nregard it as a complementary measure, while the component-wise MAE remains the most direct and\ncomprehensive evaluation metric. An additional figure sampled according to this metric is embedded in\nResponse Fig. 5 (Supplementary Fig. 1), to present a larger variety of prediction results. The\neffectiveness of this metric can be visualized directly from Response Fig. 6 (Supplementary Fig. 4),\nwhich presents the anisotropy strength as a function of photon energy, and shows good agreement\nbetween the target spectra and model predictions. Response Fig. 6(a) shows that the faithful\nreproduction of the anisotropy strength is sufficient to reconstruct the full dielectric tensor spectrum with\nhigh fidelity. Response Fig. 6(b) further illustrates this point by comparing individual tensor\ncomponents, where close agreement between the target and prediction confirms the reliability of\nMAEnorm as a compact yet informative diagnostic for evaluating anisotropic features.\naniso \nResponse Figure 5. (Supplementary Figure 1.) Representative dielectric spectra selected according to the anisotropy norm\nerror MAE0345. Targets (left) display the physically correct anisotropy, while model predictions (right) are benchmarked\n/0123 \nunder this scalar measure. Sampling by this metric provides an interpretable view of how well anisotropy strength is captured\nacross frequency and complements the component-wise MAE reported in the main text.\nTarget\nResponse Figure 6. (Supplementary Figure 4.) (a) Comparison of target (black) and predicted (red) anisotropy strength\nspectra for a representative monoclinic compound K O Sb in the test set. The close agreement indicates that recovering the\n4 8 4 \nanisotropy strength is sufficient to reproduce the full dielectric tensor spectra with high fidelity. (b) The off–diagonal 𝑥𝑧\ncomponent is nearly identical to the target, and the diagonal elements satisfy 𝑥𝑥 ≠𝑦𝑦≠𝑧𝑧, consistent with the expected\nmonoclinic symmetry. The color coding follows the same convention as in the other plots throughout this work.\n#\"!\n!\n(a) (b)\nPrediction\n#\"!\n!\nSummary of revisions made based on the preceding discussion:\n• SI, line 37—47.\n• SI Fig. 1.\n• SI, line 99—105.\n• SI Fig. 4.\nAn additional complementary measure (Question (v)): we have now included the per-component\nmetrics to report the performance under individual component more directly. The normalized MAE\n(NMAE+-) is designed to contrast the magnitude of the MAE to the characteristic magnitude, while the\nfirst derivative MAE1,+-and the KL divergence, aims to provide the shape-aware metric. These results\nare also embedded in Response Table 1 and report in the main text as table 1.\nFollowing are the details for these metrics. Again, the first metric is the mean absolute error (MAE) of\nthe dielectric tensor, together with its per–component counterpart,\n1 1\nMAE = / / 1𝜀̂ (𝜔)−𝜀 (𝜔)1, MAE+- = / 1𝜀̂ (𝜔)−𝜀 (𝜔)1.\n6𝑁 +- +- 𝑁 +- +- \n* * \n* +,- * \nBecause the magnitudes of the diagonal and off–diagonal tensor components differ significantly, we\nadditionally report a normalized error,\nMAE+-\nNMAE+-(%) = 100× ,\n𝑀+-\nwhere the normalization factor 𝑀+- is defined as\n𝑀+- = median N𝑚𝑎𝑥 𝜀 +- (𝜔)P.\n3 3 \n* \n+- \nThat is, for each material’s index 𝑖, we first take the maximum of 𝜀 (𝜔) across the frequency domain\n3 \n𝜔. We then compute the median of these maxima over all samples in the dataset. The resulting value\n𝑀+- serves as a characteristic magnitude, providing a robust measure of the typical spectral amplitude\nagainst which errors can be meaningfully normalized.\nIn addition to the absolute magnitudes, it is also important to evaluate the spectral shape of the predicted\nresponse. To this end, we compute the first-derivative MAE of each tensor component,\n𝑁 −1 \n* \n* \nwhere Δ𝜀+-(𝜔) = 𝜀+-(𝜔+𝛿𝜔)−𝜀+-(𝜔). This metric emphasizes differences in slope rather than\nmagnitude, ensuring that fine spectral features (such as peak positions and line shapes) are faithfully\nreproduced. Finally, to quantify similarity in spectral distribution from an information-theoretic\nperspective, we compute the Kullback–Leibler (KL) divergence,\n𝑝(𝜔)\n𝐷 (𝑝‖𝑞) = / 𝑝(𝜔)log ,\n45 𝑞(𝜔)\n* \nwhere 𝑝(𝜔) and 𝑞(𝜔) denote the normalized target and predicted spectra, respectively. Together, these\ncomplementary metrics evaluate not only the overall prediction accuracy but also anisotropy strength,\nspectral shape, and distributional fidelity.\nAs a complementary evaluation to address Question (iv), we also assess errors in the spherical–\nharmonic basis, where individual sub-channels map directly onto distinct anisotropy modes (e.g., 𝑥𝑥 −\n/\n𝑦𝑦 corresponds to the ℓ = 2,𝑚 = 2 channel, − (𝑥𝑥 +𝑦𝑦)+𝑧𝑧 to ℓ = 2,𝑚 = 0, and (𝑥𝑧+𝑧𝑥) to\n0 \nℓ = 2,𝑚 = 1). Detailed results for each ℓ,𝑚 sub-channels are reported in Response Table 2\n(Supplementary Table 4). We note that while these quantities offer a more symmetry-adapted view,\ntheir interpretation is less direct than in the Cartesian basis because of the mixing coefficients introduced\nduring the Cartesian–spherical transformation. This mixing slightly redistributes weight across various\nchannels, leading to small, apparent shifts relative to the Cartesian per-component MAEs.\nFrom this decomposition, two clear trends emerge. First, the isotropic channel (ℓ = 0,𝑚 = 0) exhibits\nthe smallest overall error, consistent with the dominance and smooth energy dependence of the diagonal\ntensor components. Second, the anisotropic channels (ℓ = 2,𝑚 = ±2,±1,0) show progressively larger\nMAE and NMAE values, reflecting the greater complexity and smaller absolute magnitudes of the off-\ndiagonal responses. For instance, the ℓ = 2,𝑚 = 1 channel—corresponding primarily to the 𝑥𝑧\ncomponents—shows roughly twice the error of the isotropic channel, indicating that learning these\nweaker cross-coupling terms is more challenging. Nevertheless, the predicted errors remain within a few\npercent in normalized units, confirming that the model successfully captures both isotropic and\nanisotropic tensor behavior.\nThese complementary results, together with the new visualizations and discussion provided in the SI,\nfurther demonstrate that our model not only reproduces the dominant diagonal dielectric behavior but\nalso quantitatively resolves the anisotropic and off-diagonal components that have not been previously\nlearned by non-equivariant models.\nSummary of revisions made based on the preceding discussion:\n• Main text, line 202—217.\n• Main text, Table 1.\n• SI Sec. 1., line 25—64.\n• SI Sec. 2.3., line 189—202.\n• SI Table 4.\nQuestion 4. I get the impression that the authors consider it a success that anisotropic and in particular\noff-diagonal components vanish for the more symmetric compounds. In this context, I have two\nquestions going in different directions:\n(4a) Isn't it obvious that predictions of an equivariant network preserve the symmetry.\nIn other words: Isn't it trivial that anisotropic and in particular off-diagonal components vanish for the more symmetric\ncompounds? Or am I mistaken?\nWe thank the reviewer for this insightful question. We agree that, from a physical standpoint, the\nreviewer is correct: for high-symmetry compounds, anisotropic components—particularly the off-\ndiagonal terms—are symmetry-forbidden and should vanish.\nHowever, in the context of an equivariant neural network, this outcome is not entirely trivial. Our model\nis designed to respect SO(3) rotational equivariance, meaning that if the input atomic structure is\nrigidly rotated globally, the predicted tensor rotates in the same manner. This guarantees global\ngeometric consistency but does not explicitly encode the discrete point-group symmetries of a specific\ncrystal (e.g., cubic, tetragonal). Instead, these crystal-specific symmetry relations must be learned from\nthe atomic geometry and the data distribution during training. Whether the network recovers these point-\ngroup constraints faithfully depends on two key aspects:\n(i) Numerical realization vs. exact symmetry: While equivariance ensures that the functional form of\nthe network preserves rotational consistency, the learned parameters are still determined from finite data\nand finite precision. Consequently, small but nonzero off-diagonal elements may appear early in training\nor for nearly symmetric structures. Such deviations arise from numerical noise, incomplete convergence,\nor slight asymmetries in the data rather than from a fundamental violation of symmetry. Prior works\nhave mentioned this issue and addressed this by enforcing hard symmetry masks on the outputs [a], but\nwe found that this approach limits the model’s generalization under distortions, where the effective\nsymmetry is reduced and the masks become ill-defined.\n[a] Yan, K., Saxton, A., Qian, X., Qian, X., & Ji, S. (2024). A space group symmetry informed network\nfor o (3) equivariant crystal tensor prediction. arXiv preprint arXiv:2406.12888.\n(ii) Emergent learning of crystal symmetry during convergence: As training progresses, the model\ngradually learns to recover the correct symmetry relations, which can also be directly visualized in\nprediction results. To examine this process, we applied the same shear-strain test to two unconverged\nmodels (Epochs = 1 and 10) and compared the predictions with those of the fully converged model. We\nfound that the off-diagonal components evolve systematically with training—initially underestimated by\nup to an order of magnitude, then they progressively approach the larger magnitude and functional form\n(see Response Fig. 7). This convergence behavior indicates that the model learns the underlying point-\ngroup symmetry from data, rather than inheriting it trivially from the architecture. In other words, the\ncorrect anisotropic relations emerge as a learned property of the trained model, not as an automatic\nconsequence of SO(3) equivariance alone.\nIn summary, while SO(3) equivariance provides the theoretical foundation for rotational symmetry, the\nfaithful reproduction of discrete crystal symmetries requires both exact structural symmetry in the inputs\nand proper convergence of the learned parameters. We therefore consider it important—and non-\ntrivial—that our model reproduces vanishing off-diagonal components for high-symmetry crystals and\nrecovers the expected anisotropy under structural distortions.\nEpoch = 1 Epoch = 10 Converged model\nResponse Figure 7. Evolution of the off-diagonal dielectric tensor component |𝜀2̅,.| under increasing shear angle for models\ntaken at different training epochs (Epoch = 1, 10) compared with the fully converged model. The unconverged models\n(Epoch = 1 and 10) systematically underestimate the magnitude of the response, while the fully converged model accurately\nreproduces both the amplitude and the functional dependence, closely following the sine trend expected from symmetry\nanalysis. This progressive improvement illustrates that the network learns the correct anisotropic behavior and point-group\nsymmetry relations from data, rather than inheriting them trivially from the SO(3)-equivariant architecture.\nSummary of revisions made based on preceding discussion:\n• Main text, line 268—272.\n• Main text, line 308—315.\n(4b) I like the comparison behind Supplementary Figure 1.\nWe appreciate the reviewer’s positive comment on Supplementary Fig. 1. We included this comparison\nprecisely to highlight how the model behaves across symmetry-preserving versus symmetry-breaking\ncases, and we are glad the reviewer found it helpful.\nQuestion 5. Some early papers showing the ability of GNN to predict scalar spectra have not been cited,\ne.g. M. Grunert et al, Phys. Rev. Materials 8 (2024), L122201; A. Ibrahim and C. Ataca, ACS Appl.\nMaterials & Interfaces 16 (2024), 41145-41156\nWe thank the reviewer for pointing this out. It is indeed important to acknowledge the early work\ndemonstrating the ability of GNNs to predict scalar spectra. We have now added the suggested\nreferences in the revised manuscript (cited as Refs. [27] and [31]) to properly recognize these\ncontributions in the context of our study.\nSummary of revisions made based on preceding discussion:\n• Main text, adding Refs [27], [31].\nQuestion 6. Why is section 4.2 called 'tight binding calculations'?\nI see nothing in the entire section that has anything to do with tight binding except for the general\nreferences [41,49-52] and a certain similarity between the Wannier orbitals mentioned there and the\ntight-binding picture.\nWe thank the reviewer for this point. We agree that the original title “Tight-binding calculations”\nwas potentially misleading, since the section itself focused on evaluating the dielectric function from\nfirst-principles data. Our intention was to emphasize that the procedure—extracting the Hamiltonian,\noverlap, and position operators from OpenMX and projecting them onto a localized orbital basis—is\nmathematically equivalent to constructing a Wannier-based tight-binding representation of the ab initio\nHamiltonian. This tight-binding–like representation provides the foundation for efficiently evaluating\nBrillouin-zone integrals in the Kubo formula.\nTo avoid confusion, we have revised the manuscript accordingly: the dielectric function workflow has\nbeen merged with the ab initio calculation section, which is now titled “Ab initio calculations” (Section\n4.1). In the text, we explicitly clarify that the use of localized orbitals provides a Wannier/tight-binding–\nlike framework for computing optical responses.\nSummary of revisions made based on preceding discussion:\n• Main text, line 384—416.\nQuestion 7. Many people, including https://en.wikipedia.org/wiki/Permittivity interpret 'permittivity' as\nepsilon(omega) and write epsilon = epsilon..0 * (1+chi). The present manuscript deviates from this\nusage in two ways: (i) chi is called permittivity (ii) chi includes epsilon..0 (see Eq. (13)).\nThese choices should be stated explicitly in order to avoid confusion.\nEven better would be to present the results for dimensionless quantities, such es epsilon..rel = epsilon /\nepsilon..0 or the conventionally defined chi. Most people have a feeling for dimensionless epsilon..rel\nvalues, but not for values in units of milifarad/m.\nWe thank the reviewer for this important clarification. We agree that our earlier terminology was\nnonstandard and could cause confusion. Following the reviewer’s suggestion, we have adopted the\nconventional notation:\n• 𝜒+-(𝜔) is now referred to as the (dimensionless) susceptibility.\n• 𝜀+-(𝜔) = 𝜀 ( 𝜒+-(𝜔) +𝛿+-) is the permittivity.\n6\n• Results are now reported in terms of the dimensionless relative permittivity 𝜀\n+-(𝜔)\n=\n767(*)\n,\n(\n78\n+- +- \nwhere we denote the real part as 𝜀 and the imaginary part as 𝜀 .\n/ 0 \nThis revision ensures consistency with conventional usage and allows readers to directly interpret the\nreported spectra in familiar dimensionless units.\nSummary of revisions made based on preceding discussion:\n• Main text and SI, all 𝜒 were changed to 𝜀 with 𝜀 and 𝜀 to represent the real part and imaginary\n/ 0 \npart.\nQuestion 8. The manuscript discusses only spectra calculated within the Kohn-Sham independent\nparticle picture.\nThe fact that and why these might differ dramatically from experimental spectra should at least be\nmentioned.\nWe thank the reviewer for raising this important point. Indeed, all spectra in the present work are\ncomputed within the Kohn–Sham independent-particle approximation (IPA), which neglects\nquasiparticle corrections and excitonic effects. As a result, the absolute position and detailed line shapes\nof the calculated spectra may deviate from experiment. In the revised manuscript, we have explicitly\nstated this limitation in Section 4.2.\nIt should be noted that recent work has shown that machine-learning models trained on IPA spectra can\nbe fine-tuned or transferred to RPA-level data, thereby systematically incorporating local-field effects\n[b]. This also highlights the potential advancement of our framework: while the present work is based on\nIPA data, the same model architecture could in principle be extended to higher level many-body\napproaches (e.g., GW or GW–BSE) to capture quasiparticle corrections and excitonic effects. The\nabove-mentioned reference is also added as refence [32] in the main text.\n[b] Grunert, M., Großmann, M. & Runge, E. Machine learning climbs the Jacob’s Ladder of\noptoelectronic properties. Nat Commun 16, 8142 (2025).\nSummary of revisions made based on preceding discussion:\n• Main text, line 365—373.\n• Relevant reference added in main text as [32].\nExemplary issue:\nIssue 1. The subsection \"Author contribution\" states \"T.-W. H., Z. F. and Q. Y. conceived the\nexperiments.\"\nGiven that no experiments are presented, this surely is a cut\\&paste error.\nIssue 2. The reference list is full of names of persons or elements which are not capitalized: \"weyl\nmetal\", \"van der waals\", \"from h to kr\".\nThis can easily be fixed, but should have been fixed before submission.\nIssue 3. Similarly: If the 'S' in Ref. [54] \"Sgdr: Stochastic.. \" is capitalized, one should capitalize the\nName of the Method 'SGDR' as well.\nWe thank the reviewer for pointing out these issues. The use of “experiments” in the Author\nContributions section was a cut-and-paste error, and we have corrected it to “study.” The lowercase\nreferences (e.g., weyl metal, van der waals, from h to kr) and the inconsistent capitalization of “SGDR”\nwere due to a LaTeX compilation issue in the template we used. We have carefully revised the reference\nlist and corrected these capitalization errors accordingly.\nSummary of revisions made based on preceding discussion:\n• Main text, line 490: T.-W. H., Z. F., and Q. Y. conceived the study.\nIssue 4. In the abstract, the authors state \"our model achieves a mean absolute error (MAE) of 21.181\nmillifarads per meter (mF/m)\". We teach our students that specifying an error estimate with 10^(-5)\naccuracy makes no sense at all. Since I have just quoted this sentence, I would like state again that it\nwould probably make more sense to consider dimensionless quantities. Honestly, how many readers are\nfamiliar enough with permittivity values to appreciate this numerical error value? In contrast, many\nknow that epsilon..relativ values are of order 1-10.\nIssue 5. Furthermore, in an abstract in a renowned journal such as Nature Comm., it is not necessary to\nexplain that ‘mF’ stands for ‘millifarad’ in the context of permittivity.\nWe thank the reviewer for raising these two important points. We fully agree that reporting errors in\nabsolute units such as millifarads per meter (mF/m) is potentially confusing, particularly since relative\npermittivity is by convention a dimensionless quantity. In the revised manuscript, we have corrected this\nby consistently referring to the complex relative dielectric tensor and by reporting all error metrics in\ndimensionless units.\nAccordingly, we have also simplified the abstract: the reference to “mF/m” has been removed, and the\nreported MAE values are now presented in a concise, dimensionless form, which is both physically\ncorrect and more accessible to a broad readership.\nSummary of revisions made based on preceding discussion:\n• We have revised the manuscript to consistently use “dielectric” in the main text and SI where\nrelevant.\nIssue 6. In addition to specific references to concrete points, this work naturally also contains more\ngeneral references. However, some of these seem to have been chosen extremely arbitrarily:\nFor example, the introduction begins with \"... which together form the foundation for interpreting a wide\nrange of optical phenomena in solid-state systems [1, 2].\"\nThe first reference \"[1] Lucarini, V. Kramers-Kronig relations in optical materials research (2005)\"\nrefers to a book on a very special aspect of optical spectra, without indicating the publisher.\nThe second reference \"[2] Boyd, ...Nonlinear Optics, 1097–1110 (Springer..., 2023)\" refers to a chapter\nof a book, specifying the page numbers and the name of the chapter, but not the name of the book.\nAnd why is a chapter on 'nonlinear' optics cited if the whole manuscript addresses 'linear' optics?\nI find the references for \"we selected materials ... appropriate for optoelectronic applications [34, 35]\"\nHow did ref. [34] \"Band parameters for III–V compound semiconductors and their alloys\" influence the\nchoice of the 1432 bulk semiconductors considered in this study.\nWe thank the reviewer for pointing out the inconsistencies in our references. In the revised manuscript,\nwe have standardized all book citations by including the complete title, publisher, edition, and year. For\nexample, the reference to Lucarini now cites the full book Kramers–Kronig Relations in Optical\nMaterials Research (Springer, 2005), and Boyd is cited as Nonlinear Optics, 3rd edition (Academic\nPress, 2008).\nWe also clarify why Boyd’s text is cited, despite the manuscript focusing on linear optical responses.\nSpecifically, Chapter 1 of Nonlinear Optics discusses the symmetry-based restrictions on tensor\ncomponents, including which dielectric tensor elements vanish in crystals with certain point-group\nsymmetries. These arguments are general and apply equally to linear and nonlinear optical tensors. For\nthis reason, we retained Boyd as a reference, but now explicitly cite Chapter 1 rather than an internal\npage number, to avoid confusion.\nRegarding dataset selection, we corrected the earlier reference and now clearly state that materials were\nchosen from the Materials Project [c], with band gap ranges guided by established optoelectronic\nscreening criteria [d, e].\n[c] Jain, A. et al. Commentary: The Materials Project: A Materials Genome Approach to Accelerating\nMaterials Innovation. APL Materials 1, 011002 (2013).\n[d] Kang, Y., Youn, Y., Han, S., Park, J. & Oh, C.-S. Computational Screening of Indirect-Gap\nSemiconductors for Potential Photovoltaic Absorbers. Chem. Mater. 31, 4072–4080\n[e] Yang, R. X., Horton, M. K., Munro, J. & Persson, K. A. High-throughput optical absorption spectra\nfor inorganic semiconductors (2022). Preprint at https://arxiv.org/abs/2209.02918.577.\nSummary of revisions made based on preceding discussion:\n• The relevant references in main text [3], [41], [42], and [43] have been updated accordingly.\nIssue 7. Fig. 6(b) includes as inset a linear fit of chi..xz vs. shear-strain angle.\nNaively, I would aspect a sine relation. If I look at the data, I get the feeling that a sine shape starting\nslightly deeper and then leveling off slightly could possibly fit even better. It should also be relatively\nsimple for the authors to compute the actual values of the sheared/stretched crystals, to verify whether\nthe model predicts the correct trend.\nWe thank the reviewer for this insightful comment. Indeed, for small shear-strain angles (< 10∘), the\ndata can be well approximated by a linear fit, consistent with the Taylor expansion sin𝜃 ≈ 𝜃. This is\nwhy we originally included only the linear inset fit. When the shear distortion is artificially extended to\nlarger angles (10∘ ≤ 𝜃 < 25∘), the response follows a clear sine trend, as expected. However, we note\nthat such large distortions are structurally unstable and energetically unfavorable and thus are less\nphysically relevant in actual relaxed shear-strained crystals.\nTo clarify this point, we have added the Response Fig. 8 (Supplementary Fig. 13) to include both\nregimes: (i) a linear fit in the small-angle range (0–10°), which corresponds to realistic structural\ndeformations, and (ii) a sine fit across the full range (0–25°), showing the global behavior if distortions\nare pushed beyond the physically stable regime. Both achieve high 𝑅0 values, but the sine function\ncaptures the expected asymptotic trend at large angles. We also include this effect in the discussion\nsection of the main text.\nFurthermore, we verified the model predictions against explicitly sheared crystal structures obtained\nfrom first-principles calculations. These confirm that the physically relevant regime (<10°) exhibits\nnearly linear behavior, while larger angles yield a sinusoidal dependence. Our model successfully\nrecovers both, demonstrating robustness in capturing the anisotropic response under shear strain. For\ncompleteness, we also examined uniaxial strain, again validated against first-principles data, and\nobserved the same behavior. This confirms that our model reliably captures symmetry perturbations\nacross different types of strain.\n(a) ML predicted trend DFT calculated trend\n(b)\nResponse Figure 8. (Supplementary Figure 13.) Benchmarking of model predictions against ab initio dielectric tensors\nunder symmetry-breaking perturbations. (a) Uniaxial strain (0 — 5%) applied along the 𝑐-axis of AgCl, showing a linear\ndependence of anisotropy mean along the frequency, namely, 𝛥2,,9..(𝜔) captured equally by the model and ab initio\ncalculations. (b) Shear distortions: for small shear angles (<10∘) applied on orthorhombic crystal Ag AsS , the response is\n3 4 \nnearly linear, while for extended angles (10∘ ≤𝜃 <25∘) the dependence follows a sinusoidal trend. In both cases, the model\nreproduces the expected functional form with high 𝑅; values, confirming its ability to capture anisotropic responses under\nstrain and shear perturbations.\nSummary of revisions made based on preceding discussion:\n• Main text, Fig. 6.\n• Main text, line 328–331.\n• Main text, line 343–346.\n• SI Section 5. Line 328–344.\n• SI Fig. 13.\nIssue 8. Given that an equivariant model should by definition only predict the correct tensor elements,\nthe purpose of the section on masking the resulting tensor is not clear to me.\nWe thank the reviewer for this comment. While the equivariant architecture enforces tensorial\ntransformation rules by construction, exact symmetry is not guaranteed in practice due to numerical\nfactors. The masking step was introduced to illustrate how such numerical artifacts can be identified and\nsuppressed. Further discussion of this issue is provided above under our response to Reviewer #1,\nQuestion 4.\nSummary of revisions made based on preceding discussion:\n• Main text, line 268–272.\nIssue 9. The prediction time of TSENN should also be contrasted with the time taken by the equivalent\nab initio calculation.\nWe thank the reviewer for this suggestion and acknowledge that it is important to address this point. As\na representative example, we recalculated the dielectric function of Ag AsS (Materials Project ID: mp-\n3 4 \n9538, 16 atoms per unit cell). The DFT self-consistent calculation converged in 119 s using a hybrid\nMPI/OpenMP setup on 2 AMD EPYC 7763 CPUs (256 threads in total) with the default workflow\nsettings. The subsequent dielectric tensor evaluation, carried out with Brillouin-zone integration\nparallelized over 128 processes, required 22 min 9 s. Accounting for job orchestration overhead, the full\nworkflow—from FireWorks launch to data storage and workflow termination—took 26 min 32 s.\nBy contrast, the GNN-based workflow is orders of magnitude faster. For the same material, graph\nconstruction took 6.96 ± 0.01 ms, forward inference 695.20 ± 0.01 ms, and Kramers–Kronig evaluation\n570.04 ± 0.01 ms, giving a total runtime of 1.18 s when performed on a single NVIDIA RTX 4090\nGPU. Compared to the 26 min 32 s required for the full DFT workflow, our GNN-based workflow\neffectively reduced the half-hour DFT workflow to just over a second.\nSummary of revisions made based on preceding discussion:\n• This contrasted advantage is now outline in Main text, line 456—465.\nIssue 10. Some minor questions regarding the choice of different parameters:\nWhy are the chosen batch sizes of the model so small - is the VRAM usage of the model that high?\nIf so, that should be indicated as well.\nWe thank the reviewer for this question. Indeed, batch size is one of the most sensitive hyperparameters\nin equivariant neural networks, as also reported in the literature [f]. Smaller batch sizes often yield better\ngeneralization because the model avoids over-smoothing across samples.\n[f] Batzner, S., Musaelian, A., Sun, L. et al. E(3)-equivariant graph neural networks for data-efficient\nand accurate interatomic potentials. Nat Commun 13, 2453 (2022).\nIn our study, we observed two clear regimes:\n• Too small (batch size = 1): the model sees only one data per step, which makes optimization\nunstable and reduces the statistical averaging needed for robust gradient estimates. This leads to\npoor generalization (48.32%).\n• Too large (batch size ≫ 8): the model tends to average data within a batch, washing out subtle\ndifferences between the individual materials, especially in anisotropy, where fine details are\neasily lost.\nTo quantify this, we defined 𝜌norm on the mean absolute error in anisotropy strength, together with a\naniso \nfixed threshold 𝜏 = 0.1, to evaluate how well the predicted spectra recover anisotropic behavior.\n1\n𝜌norm = e / 𝕀(MAEnorm < 𝜏)g×100,\naniso 𝑁 aniso \n;\nThe results are:\n• Batch size 1 → 48.32%\n• Batch size 2 → 62.99%\n• Batch size 4 → 70.67% (peak performance)\n• Batch size 8 → 70.46%\n• Batch size 16 → 67.25%\n• Batch size 32 → 57.26%\n• Batch size 64 → 49.16%\nThese results show that a sweet spot exists around batch sizes of 4–8, while both very small and very\nlarge batch sizes yield sub-optimal performance. For very large batch sizes (>32), the anisotropy\nbecomes less pronounced compared to smaller batches. We attribute this to the optimizer being more\nprone to local minima and averaging out sharper anisotropic signals due to stochastic smoothing.\nTherefore, smaller batch sizes are generally more favorable. To ensure the robustness of this choice, we\nrepeated the hyperparameter optimization using Optuna, explicitly scanning batch sizes from 1 to 16,\nand selected the best model based on anisotropy strength matching.\nResponse Figure 9. Variation of the anisotropy norm metric, 𝜌norm, as a function of batch size. Batch sizes from 1 to 64 were\naniso \ntested: very small batches (e.g., 1) lead to unstable training due to insufficient statistical averaging, while very large batches\ntend to oversmooth across the mini-batch, averaging out anisotropic features and losing important spectral details. A clear\noptimum is observed around batch sizes of 4–8, which balance stability and resolution. An additional Optuna scan over this\nrange confirmed that a batch size of 8 provides the best trade-off and was therefore selected for the final model.\nSummary of revisions made based on preceding discussion:\n• Main text line 444—449.\nIssue 11. Why was the broadening of the dielectric function chosen to be 50 meV, a relatively small\nvalue?\nWe thank the reviewer for this insightful question. Physically, the broadening parameter corresponds to\nthe finite relaxation time of excited carriers (𝜂 ∼ 1/𝜏), and a value of 50 meV is consistent with typical\nrelaxation times on the order of 10 femtoseconds in semiconductors [g,h,i]. Numerically, the broadening\ncontrols the resolution of the spectra obtained with a converged k-grid: smaller values yield sharper\npeaks and more pronounced spectral features. While this higher resolution makes the raw spectra more\npeaky—and therefore harder for the GNN to learn directly—subsequent smoothing through\ninterpolation or filtering mitigates this issue. Since the broadening effectively plays the same role as\nlifetime-induced spectral smearing, we expect that choosing a different but physically reasonable value\nwould not qualitatively alter our conclusions.\n[g] Del Sole, R. & Girlanda, R. Optical properties of semiconductors within the independent-particle\napproximation. Phys. Rev. B 43, 7136–7143 (1991).\n[h] N. Sivadas, S. Okamoto, and D. Xiao. Gate-controllable magneto-optic Kerr effect in layered\ncollinear antiferromagnets. Phys. Rev. Lett. 117, 267203 (2016).\n[i] Zhou, Z., Cao, G., Liu, J. et al. High-throughput prediction of the carrier relaxation time via data-\ndriven descriptor. npj Comput Mater 6, 149 (2020).\nSummary of revisions made based on preceding discussion:\n• Main text line 408—411.\n• The above-mentioned reference is now added in the main text as reference [61], [62], and [63].\nReviewer #1 (Remarks on code availability):\nWhile some code has been made openly available, the code for the ablation experiments is missing.\nWhile this is very important for possible users of TSENN, it would be beneficial for the review as well.\nFurthermore, the workflows for the generation of the optical data appears to be missing too.\nWe thank the reviewer for this important comment. In the revised submission, we have now made the\ncode for the ablation experiments openly available on the GitHub, together with detailed checkpoints to\nensure reproducibility. We have also added the workflow used for generating the optical data, which is\navailable at: https://github.com/qmatyanlab/OpenMX-workflow. This repository contains the full\npipeline for data preparation and optical response generation from OpenMX, allowing the entire\nworkflow to be independently reproduced and extended by the community.\nSummary of revisions made based on preceding discussion:\nMain text, line 466: Data availability\nAll datasets and optimized model weights used in this study are openly available at\nhttps://github.com/qmatyanlab/TSENN. The repository also includes reproducible Jupyter notebooks\ndemonstrating data usage and reproducing the figures presented in this work.\nMaine text, line 470: Code availability\nThe third-party code OpenMX used in this study is publicly available at https://www.openmx-\nsquare.org. The custom workflows developed to compute the dielectric functions are accessible at\nhttps://github.com/qmatyanlab/OpenMX-workflow and were implemented using the Julia programming\nlanguage https://julialang.org.\nReviewer #3 (Remarks to the Author):\nHsu et al. have demonstrated a powerful neural network architecture, namely TSENN, that aims to relate\nthe crystal structure to the full anisotropic frequency-dependent dielectric permittivity. The demonstrated\narchitecture shows high prediction accuracy, confirmed by the KK relations with capabilities of\nmaintaining the anisotropy. Furthermore, authors have presented two case studies of strain-induced\nsymmetry breaking and the resulting change in dielectric permittivity.\nThe data-driven results are timely and interesting, and particularly useful for materials screening for\nboth new materials discovery as well as application purposes – a helpful methodology for both academia\nand industry. However, the advantage of simplifying the symmetry processing and accuracy in\npermittivities, particularly below/near bandgap, would require further clarification.\nAuthor response to Reviewer #3:\nWe thank the reviewer for their thoughtful and positive evaluation of our work and for recognizing the\ncontributions of TSENN to the field. We appreciate the reviewer’s comments regarding the model\naccuracy near and below the bandgap. In response, we have carefully examined and discussed the\nprediction results in this regime, specifically analyzing the spectra around the plasma frequency and the\noptical bandgap extracted from Tauc plots. These analyses confirm that the model remains physically\nconsistent and accurately reproduces dielectric trends across the bandgap region. The corresponding\ndiscussion and figures have been added to the revised manuscript to address this point in detail as\ndiscussed below.\nDetailed responses to reviewer questions are given below:\nQuestion 1. What are the advantages of using spherical harmonic decomposition compared with using\ndielectric permittivity? For non-magnetic materials, both approaches yield the same number of variables\nfor all seven lattices, given symmetric off-diagonal dielectric components (i.e., 𝑥 = 𝑥 ) and rotational\n3< 3< \nsymmetry. Thus, to me, spherical harmonic decomposition does not offer significant differences in the\nanisotropy results and improvement in reducing the computational complexity. For example, a uniaxial\nsystem has 2 unique dielectric permittivities (xx, zz), 4 for monoclinic (xx, xz, yy, zz), and 6 for\ntriclinic, and the same variables are generated from spherical harmonic decomposition. If no significant\nadvantages, what’s new regarding this approach?\nWe thank the reviewer for this insightful question. The main advantage of spherical harmonic\ndecomposition is not in reducing the number of variables—indeed, the count of independent components\nis the same—but in providing a physics-inspired, symmetry-adapted representation. Cartesian\ncomponents do not intrinsically guarantee rotational symmetry: see Response Fig. 10 (Supplementary\nFig. 7). Direct training in the Cartesian basis leads to broken isotropy even for cubic materials, where\n𝑥𝑥 = 𝑦𝑦 = 𝑧𝑧 should hold. In contrast, the spherical representation enforces these fundamental relations\nby construction, embedding physical constraints directly into the representation without extra\ncomputational overhead.\nA second advantage is robustness under distortions. In practice, real materials rarely preserve ideal\nsymmetry—structural relaxations, externally applied strain, or thermal fluctuations can all break lattice\nsymmetry. In the Cartesian basis, tensor components mix under rotations in a complicated, non-\ntransparent way. The spherical harmonic decomposition, however, handles this mixing naturally and\nequivariantly. This is why we explicitly validated our model under strained lattices, where it maintained\nconsistent performance.\nIn summary, the spherical harmonic decomposition offers:\n(i) a symmetry-adapted representation that enforces correct physical constraints;\n(ii) natural generalization to distorted structures through SO(3)-equivariant mixing; and\n(iii) a physically interpretable framework that cleanly separates isotropic (ℓ = 0) and anisotropic (ℓ =\n2) contributions.\nThese aspects go beyond what Cartesian channels alone can provide and are central to the novelty of our\napproach.\nResponse Figure 10. (Supplementary Figure 7.) Comparison of target (left panels) and predicted (right panels) imaginary\ndielectric tensor spectra for selected cubic crystals. In the target data, isotropy enforces 𝜀̂,, =𝜀̂-- =𝜀̂.,, resulting in a single\n; ; ; \nline for the diagonal components (color-coded curves overlap perfectly). By contrast, the predictions visibly separate into\nthree distinct curves, indicating that isotropy is not preserved and the equivariance constraint is broken. This discrepancy\nhighlights the necessity of a spherical-harmonic representation to enforce rotational symmetry and maintain consistency with\ncubic lattice constraints.\nSummary of revisions made based on preceding discussion:\n• SI Section 2.4, line 203—232.\n• SI Fig. 7.\nQuestion 2. Typically, the dielectric permittivity below the plasma frequency matters more as that\nhighlights distinct material characteristics. Looking at figure 4b, the spectrum shows more deviation in\nthe absorption spectrum, while better agreement above plasma frequency (where fewer material\ndistinctions will be shown).\na. Taking more frequencies above the plasma frequency reduces the statistical weight of the absorption\nspectrum; thus, I think authors should have a statistical discussion and figure focusing on frequencies\nbelow the plasma frequency.\nWe thank the reviewer for this helpful suggestion. In the revised manuscript we benchmark the dielectric\ntensor separately below and above the plasma frequency 𝜔 (extracted from the energy-loss spectra) and\n= \nreport both MAE and NMAE per component.\nAs shown in the Response Fig. 11 (Supplementary Fig. 12), errors are markedly larger below 𝜔 ,\n= \nwhere material-specific absorption features dominate:\n• MAE (below 𝜔 ): diagonals ~ 0.489–0.532; off-diagonals ~ 0.146—0.183.\n= \n• MAE (above 𝜔 ): diagonals ~ 0.043–0.050; off-diagonals ~ 0.012—0.015.\n= \n• NMAE (below 𝜔 ): components range ~ 6%–13%.\n= \n• NMAE (above 𝜔 ): components drop to ~ 0.6%–1.1%.\n= \nThis frequency-resolved analysis confirms the physical picture: below 𝜔 the spectra carry strong,\n= \nmaterial-dependent oscillator strengths and are harder to predict, while above 𝜔 the response is\n= \nsmoother and less discriminative, yielding much smaller errors.\nImportantly, compared to the aggregate metrics in Response Table 1, where all frequencies are pooled,\nconditioning on 𝜔 reveals nearly a two time increase in the measured errors for the low-frequency\n= \n(absorption) region. This avoids the dilution effect from having more frequency points at high energies\nand directly addresses the reviewer’s concern about statistical weighting.\nWe considered training with an energy-resolved mask to up-weight 𝜔 < 𝜔 , but this would require prior\n= \nknowledge of 𝜔 for each material and could introduce material-dependent biases at inference time.\n= \nSince 𝜔 varies widely across compounds, we instead adopted an energy-independent training loss and\n= \nreport frequency-resolved metrics, which we believe provide a cleaner and more general solution. The\nnew figure and discussion have been added accordingly.\n(a)\n(b)\nResponse Figure 11. (Supplementary Figure 12.) Per-component MAE*+ and NMAE*+ of the imaginary dielectric tensor,\nevaluated separately below and above the plasma frequency 𝜔 . Once the spectra are partitioned at 𝜔 , the errors become\n< < \nnearly twice as large compared to the aggregate values reported in Table 1 of the main text. The low-frequency regime (𝜔<\n𝜔 ) dominates the error, with significant contributions from both diagonal and off-diagonal components, whereas the high-\n< \nfrequency regime (𝜔>𝜔 ) exhibits substantially smaller errors, with NMAE*+ reduced to only a few percent per\n< \ncomponent. This trend is consistent with the physical picture: below 𝜔 , the dielectric response captures rich, material-\n< \nspecific absorption features with pronounced oscillator strengths, making this regime intrinsically more sensitive to\npredictive errors, while the smoother high-frequency response is easier to approximate.\nSummary of revisions made based on preceding discussion:\n• Main text, line 354—364.\n• SI, line 309—327.\n• SI Fig. 12.\nb. Accurately predicting the bandgap has been a long-standing topic in the community. It would be\nbeneficial if authors could compare the bandgap from predicted dielectric permittivities (e.g., Tauc\nmethod[j]) vs the calculated bandgap to highlight the prediction accuracy.\n[j] B. D. Viezbicke, S. Patel, B. E. Davis, D. P. Birnie III, Phys. Status Solidi B 2015, 252, 1700.\nWe thank the reviewer for this helpful suggestion. To directly assess the accuracy of our bandgap\npredictions, we extracted optical bandgaps from both the DFT-calculated and ML-predicted dielectric\nfunction using the Tauc method [j]. By definition, the optical bandgap corresponds to the lowest photon\nenergy that can be absorbed to produce an allowed electronic excitation. Consequently, if direct optical\ntransitions are forbidden at the electronic band edge, the extracted optical bandgap will necessarily be\nlarger than the fundamental electronic bandgap. This explains why in practice optical bandgaps are\nobserved to be greater than or equal to the electronic bandgap, making the comparison between the two\nphysically consistent.\nThe left panel of the Response Fig 12. (Supplementary Fig. 11) compares the DFT electronic\nbandgaps with optical bandgaps obtained from Tauc fits to the calculated dielectric functions, showing\nthe expected trend. The right panel compares Tauc-derived optical bandgaps from the ML-predicted\ndielectric spectra against the reference Tauc-derived values, achieving an 𝑅0 of 0.897 on the test set.\nThis demonstrates that our model not only reproduces imaginary dielectric spectra with high fidelity but\nalso encodes sufficient information to recover the full dielectric function, and consequently the\nabsorption spectra and fundamental bandgaps, with high accuracy.\n(a) (b)\nResponse Figure 12. (Supplementary Figure 11.) (a) Comparison of DFT-calculated electronic bandgaps with optical\nbandgaps obtained from Tauc fits. (b) Comparison of target Tauc-derived optical bandgaps with values extracted from ML-\npredicted absorption spectra, showing strong agreement (𝑅; =0.897).\nSummary of revisions made based on preceding discussion:\n• SI, line 295—308.\n• SI Fig. 11.\nQuestion 3. Magnetic materials can have non-symmetrical dielectric permittivity. Would current\nmethods capture all 9 dielectric permittivities without assuming symmetrical components?\nWe thank reviewer for this insightful question. The current method can be extended to capture all 9\nindependent components of the dielectric tensor given the appropriate dataset. In the present work we\nenforced the usual 𝑖𝑗 = 𝑗𝑖 symmetry valid for non-magnetic crystals, which reduces the output to the six\nunique upper-triangular components. For magnetic or magneto-optical materials, one can simply lift this\nconstraint and predict the full 3×3 tensor. The additional three antisymmetric channels then naturally\ndescribe the magneto-optical contributions (e.g., 𝜀!> ≠ 𝜀>!), without further modification of the\nframework. In fact, these antisymmetric terms can be mapped directly onto spherical tensor components,\nfor example:\n1 1 1\n𝑌 = (𝜀\"! −𝜀!\" ),𝑌 = (𝜀!> −𝜀>! ),𝑌 = (𝜀>\" −𝜀\">),\n/,?/ /,6 /,/\n√2 √2 √2\nwhich makes explicit how the antisymmetric part of the tensor enters in the magneto-optical response. In\nthe e3nn notation, this corresponding to adding an additional channel, namely 𝑁 0𝑒+𝑁 1𝑜+𝑁 2𝑒,\n* * * \nat each frequency 𝑁 , you will have 0𝑒+1𝑜+2𝑒 = 1+3+5 total 9 components.\n* \nSummary of revisions made based on preceding discussion:\n• SI, line 428—436.\nReviewer #4 (Remarks to the Author):\nThis paper presents a model that predicts the full, frequency-dependent dielectric tensor directly from\nstructure using an equivariant GNN (E(3)NN). The outputs are decomposed into ℓ=0 (isotropic) and ℓ=2\n(anisotropic) spherical components. The dataset includes 1,432 non-magnetic semiconductors from the\nMaterials Project, recalculated with OpenMX (numerical pseudo-atomic orbitals). The imaginary part of\nthe dielectric tensor is computed via the Kubo–Greenwood formalism for 0–30 eV (0.01 eV step, but\ndownsampled by interpolating the data onto a 300-point photon energy grid for GNN training). The\nauthors report a test MAE of 21.181 mF/m and use Kramers–Kronig to get the real part. They also show\nsymmetry checks under strain/shear that indicate the model respects crystallographic constraints. I\nsuggest publication of this manuscript after minor revision.\nAuthor Response to Reviewer #4:\nWe sincerely thank the reviewer for the positive evaluation of our manuscript and for recommending\npublication after minor revisions. We are grateful for the constructive comments, which have helped us\nfurther sharpen the novelty of our work, strengthen the validation, and provide additional context.\nBelow we respond to each point in detail.\nDetailed responses to reviewer questions are given below:\nQuestion 1. In the intro, the authors write:\n“Notably, existing studies mostly focus exclusively on the dielectric constant, i.e., the ω → 0 static limit.\nIn contrast, the joint prediction of tensorial structure and spectral behavior—specifically, the full\nfrequency-dependent dielectric tensor—remains a largely unexplored area. In this work, we address this\ngap by leveraging the expressive power of higher-order components within the equivariant GNNs to\ndirectly predict the full tensorial, frequency-dependent permittivity, enabling a more comprehensive\nunderstanding of the anisotropic optical response of crystalline materials.”\nI think this part would benefit from a quick mention of the literature. It might be fair to say that for\ntensorial spectra, the area has been explored before but more indirectly and in a fundamental way, for\nexample in equivariant machine-learning Hamiltonians (https://www.nature.com/articles/s41524-023-\n01130-4, https://www.nature.com/articles/s41467-023-38468-8), but not specifically for the frequency-\ndependent dielectric spectrum. There *are* already some papers that address crystal structure-to-\nspectrum predictions for the frequency-dependent dielectric function, absorption spectra, or solar\nabsorbance using equivariant or invariant GNNs\n(https://pubs.acs.org/doi/full/10.1021/acsami.4c07328, https://advanced.onlinelibrary.wiley.com/doi/full/\n10.1002/adma.202409175). I think these deserve a brief mention up front, both to acknowledge the\nprevious work and to clarify that your contribution is the full tensorial extension using an equivariant\nGNN, just to sharpen the novelty claim.\nWe thank the reviewer for pointing this out and fully agree that prior literature on structure-to-spectrum\nprediction and equivariant Hamiltonians deserves acknowledgment. In the revised introduction, we have\nincorporated these works to provide a clearer context and sharpen the novelty claim of our study.\nSpecifically, we now mention GNNOpt, Optimate, and related models that predict scalar spectral\nproperties (e.g., dielectric functions, absorption, solar absorbance), as well as AnisoNet and DTNet for\nstatic tensorial properties. We also explicitly note recent advances in equivariant machine-learning\nHamiltonians, such as HamGNN, DeepH-E3, which capture tensorial spectra indirectly but still require\nexplicit Brillouin-zone integration. This makes clear that, while these works laid important groundwork,\nthe direct prediction of the full frequency-dependent dielectric tensor using an equivariant GNN—as\npresented here—remains unexplored and constitutes the key contribution of our work.\nSummary of revisions made based on preceding discussion:\n• Main text, Line 52—58.\n• Main text, Line 66—67.\n• The reference mentioned is now embedded as [25, 26, 27, 30] in the main text.\nQuestion 2. MAE alone as an error metric can miss important spectral details. Two predictions might\nhave the same MAE, but differ a lot in terms of peak positions, heights, band gap onset, or static-limit\nvalues. Please consider adding one or two spectrum shape-aware metrics on the test set (for example,\nfirst-derivative MAE, or a shape distance like earth-mover distance or KL divergence on a PDF-\nnormalized spectrum, or other relevant metrics) and report these in a small table or figure. This would\nback up validation about spectral fidelity.\nWe thank reviewer for this insightful question. We have added the MAE, per component MAE,\nNormalized MAE and two shape-aware metrics, the first-derivative MAE and KL divergence on a PDF-\nnormalized spectrum. However, we would like to point out the KL divergence cannot be applied on the\noff-diagonal component, as the value can be either positive or negative, which makes the KL ill defined.\nThe corresponding table has been embedded together with Reviewer #1 as Response Table 1. Also, the\ntable is embedded in the main text as Table 1.\nSummary of revisions made based on preceding discussion:\n• Main text, Table 1.\n• Main text, line 219—243.\n• SI Sec. 1., line 25—64.\n• SI Fig. 1.\nQuestion 3. Since the outputs are split into ℓ=0/ℓ=2, it might be helpful to see/analyze an ablation where\nyou train only on ℓ=0 and then only on ℓ=2. This may show how much accuracy comes from the\nanisotropy vs. isotropy channels, and whether there’s any leakage between channels.\nWe thank the reviewer for this excellent suggestion. As shown in Response Table 3 (Supplementary\nTable 3), training only on one channel confirms that there is no leakage: when the model is trained with\nℓ=0 only, the ℓ=2 contribution vanishes, and vice versa. This demonstrates that the isotropic and\nanisotropic channels are cleanly separated in our framework.\nResponse Table 3 reports the mean absolute error calculated in the spherical-harmonics basis, which\nprovides a natural separation into the invariant trace (ℓ = 0) and the traceless anisotropic part (ℓ = 2) of\nthe dielectric tensor. Namely,\n𝑁 B B\n* \nB * \n𝑁 B B\n* \n* \nwhere 𝑚 runs over the five components of the ℓ = 2 subspace. The total MAE can then be written as\n0\nso that the channel-wise contributions are given by\n6\n0\nBy construction,\nMAEtotal = Contrib ℓA6 + Contrib ℓA0,\nwhich allows us to quantify the relative error contributions of the isotropic and anisotropic channels\nwhile ensuring they sum exactly to the total.\nIn terms of accuracy, the ablation highlights that both channels contribute comparably to the total MAE\n(approximately 52.2–58.9% from ℓ=0 and 41.1–47.8% from ℓ=2, depending on cubic vs. non-cubic\nsystems). For cubic crystals, as expected, the error comes exclusively from the ℓ=0 channel, consistent\nwith symmetry considerations. These results confirm that the isotropic and anisotropic channels are both\nnecessary for capturing general optical responses, while remaining independent without cross-channel\nleakage.\nSummary of revisions made based on preceding discussion:\n• SI, Sec. 2.3, line 160—188.\n• SI Table 3.\nReviewer #4 (Remarks to the code availability):\nMy comments are reconciled and given in \"comments for authors\"\nWe thank the reviewer for the constructive feedback and for confirming the reconciliation of the code\navailability remarks.\nReviewer #5 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of\nthe Nature Communications initiative to facilitate training in peer review and to provide appropriate\nrecognition for Early Career Researchers who co-review manuscripts.\nWe sincerely thank the reviewer for their thoughtful input and for contributing through the Nature\nCommunications co-review program for Early Career Researchers.\nReviewer #5 (Remarks to the code availability):\nThe TSENN code is accompanied by a clear README, includes necessary scripts and notebooks, and\nappears to be a well-organized, reproducible resource for the community, although I have not personally\nrun it.\nWe thank the reviewer for their positive assessment of the TSENN code and accompanying\ndocumentation.\nSincerely,\nQimin Yan\non behalf of all co-authors\nDepartment of Physics and\nQuantum Materials and Sensing Institute\nNortheastern University\nLetter to Reviewers\nDear Reviewers,\nWe would like to begin by sincerely thanking you for the time and effort you have invested in\nreviewing and handling our manuscript. We greatly appreciate the detailed and knowledgeable\nreports, which raised important points ranging from technical clarifications in figures and\ncaptions to broader suggestions for strengthening the framing and future directions of this\nresearch. Your constructive feedback has been invaluable in helping us substantially improve the\nclarity, rigor, and overall presentation of the manuscript.\nIn the following, we address all comments in point-by-point detail: reviewer comments are\nreproduced in black text, our detailed responses are provided in red text, and the corresponding\nrevisions are marked in blue, with line number referring to the compiled PDF and highlighted in\nred within the manuscript for easy reference. In addition to implementing the requested changes,\nwe also took the opportunity to correct minor typographical errors and to standardize references\nto further improve the readability and consistency of the manuscript.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nMy young co-reviewer and I appreciate the many changes made to the manuscript based of our\nsuggestions and those of the other reviewers.\nWhile most technical shortcomings have been addressed, some remain. Thus, the decision to\naccept this manuscript for publication in Nature Communications or transfer it to a more\nspecialized journal is primarily one of editorial policy.\nWe thank the reviewer for their assessment and appreciate their overall positive evaluation of the\nrevised manuscript.\nFig. 5 best illustrates why I find it difficult to recommend acceptance: obviously, any \"prediction\nof tensorial properties\" should predict off-diagonal tensor components. These are present only in\nthe upper two rows. If the predicted (small) off-diagonal components were replaced by its\nopposite [\\epsilon -> (-1)*\\epsilon] the similarity to the target would not be much\nworse, as a simple glance at the figure shows.\nMy interpretation is that the authors do the right thing and implement it correctly, but that there\nis simply not enough training data with strong deviations from the isotropic, orientation-averaged\nresponse. Thus, the latter including the general frequency-dependent amplitude envelope has\nbeen learned, but not the former.\nWe thank the reviewer for this insightful observation. The reviewer is correct that the off-\ndiagonal dielectric tensor components in our dataset are typically small for most semiconductors,\nreflecting the near-isotropic optical response of many Materials Project structures. As a result,\nthe model learns the dominant isotropic amplitudes with higher accuracy, while the\ncomparatively rarer anisotropic patterns are more challenging to capture. We agree that this\nbehavior reflects a limitation of the nature of our dataset rather than of the model architecture\nitself.\nTo clarify this point for readers, we have added a brief discussion in the manuscript noting that\nthe distribution of tensor components in the training data influences the relative learning\naccuracy of isotropic versus anisotropic features. We believe this addition improves transparency\nregarding the scope and limitations of the current dataset, without altering the main conclusions\nof the work. Additionally, we have included an auxiliary dataset, as the primary dataset contains\nfew triclinic and monoclinic structures to increase their representation. Specifically, we restrict\nthe crystal system exclusive to triclinic and monoclinic structures and expand the filters to a band\ngap of 0.3-5 eV, an energy-above-hull range of 0-0.1 eV, compositions with 1-3 elements, and\nnumber of atomic sites from 1 to 40, yielding a total of 1,658 materials, consisting of 1,177\nmonoclinic systems and 481 triclinic systems. The larger energy-above-hull window and\nallowance for bigger unit cells are key: many low-symmetry materials are metastable or\ninherently have larger cells and would otherwise be excluded under the stricter criteria of the\nmain dataset. These relaxed constraints improve the coverage of low-symmetry systems but alter\nthe structural and chemical distributions. For this reason, the auxiliary dataset is evaluated\nseparately from the original benchmark, as it contains more metastable structures and differs\nsystematically in symmetry, complexity, and size.\nResponse Table 1 (embedded as Supplementary Table 6) summarizes the performance of our\nmodel on the auxiliary dataset. Relative to the main dataset results shown in Response Table 2\n(Table 1 of the main text), the auxiliary dataset achieves comparable overall performance but\nexhibits substantial improvements in the off-diagonal metrics, including the reduced 𝑀𝐴𝐸!\" and\n𝑁𝑀𝐴𝐸!\"(%). We further visualize the per-structure prediction accuracy stratified by the\n𝑀𝐴𝐸 distribution in Response Figure 1 (Supplementary Figure 16) and by the 𝑀𝐴𝐸$'()\n#$%&'\ndistributions in Response Figure 2 (Supplementary Figure 17). The off-diagonal components\nare predicted significantly more accurately than in the main dataset. This improvement is\nexpected, as the auxiliary benchmark effectively acts as a data augmentation step—training\nexclusively on materials that exhibit non-zero anisotropy. These results demonstrate that, when\nprovided with an appropriately anisotropy-rich dataset, our model can reliably learn the\nunderlying symmetry relationships encoded in the crystal structure and map them to the\ncorresponding tensorial spectra. The auxiliary dataset and the associated model checkpoints have\nalso been made publicly available in our GitHub repository.\nResponse Figure 1 (Supplementary Figure 16). Distribution of full-tensor MAE values across the test\nset of the auxiliary dataset (left) and representative spectra obtained by stratified sampling across quartiles\n(right). From the cumulative kernel-density-estimator (KDE) plot, we randomly selected two systems in\neach quartile (𝑄 −𝑄 ), resulting in 8 representative examples that span the full error distribution.\n! \"\nResponse Figure 2 (Supplementary Figure 17). Distribution of full-tensor MAE values across the test\nset of the auxiliary dataset (left) and representative spectra obtained by stratified sampling across quartiles\n(right). From the cumulative kernel-density-estimator (KDE) plot, we randomly selected two systems in\neach quartile (𝑄 −𝑄 ), resulting in 8 representative examples that span the full error distribution.\n! \"\nSummary of revisions made based on preceding discussion:\n• Main text, page 13, line 352—358.\n• Supplementary Section 8. Line 411 – 438.\n• Supplementary Table 6.\n• Supplementary Figure 16.\n• Supplementary Figure 17.\nThe whole paper and Supplemental Information are quite technical with a lot of definitions\nrepeated for convenience of the reader from original sources.\nWe thank the reviewer for the observation. Given the technical nature of the work, we\nintentionally included brief definitions from standard sources to keep the exposition self-\ncontained and accessible to a broad readership. We hope this improves clarity for the broad\nreadership of the journal readers who may not be familiar with some aspects of the underlying\nformalism.\nIf am not mistaken, the statement \"The dielectric tensor can be expressed as the outer product of\ntwo rank-1 vectors\" on page 3, rephrased as Eq. (1), is simply incorrect. Any tensor can trivially\nbe written as a sum of outer products.\nFurthermore, discussing these outer products is superfluous because it only leads to equations (3)\nthrough (5), which describe the action of the e3nn package, i.e., Ref. [21]. I strongly recommend\nstreamlining or fixing this part before publishing the manuscript.\nIf I misunderstood something, the authors should try to prevent that others make the same\nmistake.\nWe thank the reviewer for pointing out the ambiguity in our original phrasing. We agree that a\ngeneral dielectric tensor cannot be expressed as a single outer product of two vectors, and that\nour wording could be misinterpreted as implying a rank-1 structure. In response, we have\nremoved this part of the main text, as it was superfluous, and now refer the reader to\nSupplementary Section 6 for a detailed derivation. We have revised the text accordingly to\neliminate this ambiguity and to streamline the connection to the spherical harmonic\nrepresentation, thereby shortening this part of the discussion.\nSummary of revisions made based on preceding discussion:\n• Main text, page 3, line 105— page 4, line 110.\nA very minor suggestions is to define \\sigma_1,2 in Eq. (7) with a sentence like “where \\sigma_1\nand \\sigma 2 are learnable parameters….” [Taken from page 7 of the SI].\nWe thank the reviewer for the helpful suggestion. We have added the corresponding changes for\n𝜎 and 𝜎 .\n* +\nSummary of revisions made based on preceding discussion:\n• Main text, page 6, line 173.\nIncidentally, the respective equation on page 7 of the SI is unnumbered, unlike almost all other\nequations. On purpose?\nRegarding the Supplemental Information, we thank the reviewer for noting the inconsistency in\nequation numbering; this was unintentional, and we have corrected it in the revised SI.\nSummary of revisions made based on preceding discussion:\n• Supplementary Information, page 7, line 142.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nWe sincerely thank Reviewer #2 for their thoughtful evaluation and for participating in the\nNature Communications co-review initiative. We appreciate the time and care invested in\nassessing our work.\nReviewer #3 (Remarks to the Author):\nThe authors have addressed my concerns, and thus I support the publication.\nWe thank Reviewer #3 for their positive feedback and for acknowledging that our revisions have\naddressed the earlier concerns. We appreciate the support for publication.\nReviewer #4 (Remarks to the Author):\nWe thank the authors taking their time to revise the manuscript according to reviewers'\ncomments. We believe that the manuscript deserves publication in its current form.\nWe sincerely thank Reviewer #4 for their careful reading and for recognizing the improvements\nmade in our revision. We appreciate the supportive assessment of the manuscript in its current\nform.\nReviewer #5 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nWe thank Reviewer #5 for their contribution through the Nature Communications co-review\ninitiative and for taking the time to evaluate our revised manuscript. We greatly appreciate the\nfeedback provided.\nSincerely,\nQimin Yan\non behalf of all co-authors\nDepartment of Physics and\nQuantum Materials and Sensing Institute\nNortheastern University"}, {"round": 2, "reviewer_comments": {"reviewer_1": "My young co-reviewer and I appreciate the many changes made to the manuscript based of our suggestions and those of the other reviewers.\nWhile most technical shortcomings have been addressed, some remain. Thus, the decision to accept this manuscript for publication in Nature Communications or transfer it to a more specialized journal is primarily one of editorial policy.\nFig. 5 best illustrates why I find it difficult to recommend acceptance: obviously, any \"prediction of tensorial properties\" should predict off-diagonal tensor components. These are present only in the upper two rows. If the predicted (small) off-diagonal components were replaced by its opposite [\\epsilon -> (-1)*\\epsilon] the similarity to the target would not be much worse, as a simple glance at the figure shows.\nMy interpretation is that the authors do the right thing and implement it correctly, but that there is simply not enough training data with strong deviations from the isotropic, orientation-averaged response. Thus, the latter including the general frequency-dependent amplitude envelope has been learned, but not the former.\nThe whole paper and Supplemental Information are quite technical with a lot of definitions repeated for convenience of the reader from original sources.\nIf am not mistaken, the statement \"The dielectric tensor can be expressed as the outer product of two rank-1 vectors\" on page 3, rephrased as Eq. (1), is simply incorrect. Any tensor can trivially be written as a sum of outer products.\nFurthermore, discussing these outer products is superfluous because it only leads to equations (3) through (5), which describe the action of the e3nn package, i.e., Ref. [21]. I strongly recommend streamlining or fixing this part before publishing the manuscript.\nIf I misunderstood something, the authors should try to prevent that othersmake the same mistake.\nA very minor suggestions is to define \\sigma_1,2 in Eq. (7) with a sentence like “where \\sigma_1 and \\sigma 2 are learnable parameters….” [Taken from page 7 of the SI].\nIncidentally, the respective equation on page 7 of the SI is unnumbered, unlike almost all other equations. On purpose?\n(Remarks on code availability)", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.\n(Remarks on code availability)", "reviewer_3": "The authors have addressed my concerns, and thus I support the publication.\n(Remarks on code availability)", "reviewer_4": "We thank the authors taking their time to revise the manuscript according to reviewers' comments. We believe that the\nmanuscript deserves publication in its current form.\n(Remarks on code availability)", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.\n(Remarks on code availability)"}, "author_response": "Letter to Reviewers\nDear Reviewers,\nWe would like to begin by sincerely thanking you for the time and effort you have invested in\nreviewing and handling our manuscript. We greatly appreciate the detailed and knowledgeable\nreports, which raised important points ranging from technical clarifications in figures and\ncaptions to broader suggestions for strengthening the framing and future directions of this\nresearch. Your constructive feedback has been invaluable in helping us substantially improve the\nclarity, rigor, and overall presentation of the manuscript.\nIn the following, we address all comments in point-by-point detail: reviewer comments are\nreproduced in black text, our detailed responses are provided in red text, and the corresponding\nrevisions are marked in blue, with line number referring to the compiled PDF and highlighted in\nred within the manuscript for easy reference. In addition to implementing the requested changes,\nwe also took the opportunity to correct minor typographical errors and to standardize references\nto further improve the readability and consistency of the manuscript.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nMy young co-reviewer and I appreciate the many changes made to the manuscript based of our\nsuggestions and those of the other reviewers.\nWhile most technical shortcomings have been addressed, some remain. Thus, the decision to\naccept this manuscript for publication in Nature Communications or transfer it to a more\nspecialized journal is primarily one of editorial policy.\nWe thank the reviewer for their assessment and appreciate their overall positive evaluation of the\nrevised manuscript.\nFig. 5 best illustrates why I find it difficult to recommend acceptance: obviously, any \"prediction\nof tensorial properties\" should predict off-diagonal tensor components. These are present only in\nthe upper two rows. If the predicted (small) off-diagonal components were replaced by its\nopposite [\\epsilon -> (-1)*\\epsilon] the similarity to the target would not be much\nworse, as a simple glance at the figure shows.\nMy interpretation is that the authors do the right thing and implement it correctly, but that there\nis simply not enough training data with strong deviations from the isotropic, orientation-averaged\nresponse. Thus, the latter including the general frequency-dependent amplitude envelope has\nbeen learned, but not the former.\nWe thank the reviewer for this insightful observation. The reviewer is correct that the off-\ndiagonal dielectric tensor components in our dataset are typically small for most semiconductors,\nreflecting the near-isotropic optical response of many Materials Project structures. As a result,\nthe model learns the dominant isotropic amplitudes with higher accuracy, while the\ncomparatively rarer anisotropic patterns are more challenging to capture. We agree that this\nbehavior reflects a limitation of the nature of our dataset rather than of the model architecture\nitself.\nTo clarify this point for readers, we have added a brief discussion in the manuscript noting that\nthe distribution of tensor components in the training data influences the relative learning\naccuracy of isotropic versus anisotropic features. We believe this addition improves transparency\nregarding the scope and limitations of the current dataset, without altering the main conclusions\nof the work. Additionally, we have included an auxiliary dataset, as the primary dataset contains\nfew triclinic and monoclinic structures to increase their representation. Specifically, we restrict\nthe crystal system exclusive to triclinic and monoclinic structures and expand the filters to a band\ngap of 0.3-5 eV, an energy-above-hull range of 0-0.1 eV, compositions with 1-3 elements, and\nnumber of atomic sites from 1 to 40, yielding a total of 1,658 materials, consisting of 1,177\nmonoclinic systems and 481 triclinic systems. The larger energy-above-hull window and\nallowance for bigger unit cells are key: many low-symmetry materials are metastable or\ninherently have larger cells and would otherwise be excluded under the stricter criteria of the\nmain dataset. These relaxed constraints improve the coverage of low-symmetry systems but alter\nthe structural and chemical distributions. For this reason, the auxiliary dataset is evaluated\nseparately from the original benchmark, as it contains more metastable structures and differs\nsystematically in symmetry, complexity, and size.\nResponse Table 1 (embedded as Supplementary Table 6) summarizes the performance of our\nmodel on the auxiliary dataset. Relative to the main dataset results shown in Response Table 2\n(Table 1 of the main text), the auxiliary dataset achieves comparable overall performance but\nexhibits substantial improvements in the off-diagonal metrics, including the reduced 𝑀𝐴𝐸!\" and\n𝑁𝑀𝐴𝐸!\"(%). We further visualize the per-structure prediction accuracy stratified by the\n𝑀𝐴𝐸 distribution in Response Figure 1 (Supplementary Figure 16) and by the 𝑀𝐴𝐸$'()\n#$%&'\ndistributions in Response Figure 2 (Supplementary Figure 17). The off-diagonal components\nare predicted significantly more accurately than in the main dataset. This improvement is\nexpected, as the auxiliary benchmark effectively acts as a data augmentation step—training\nexclusively on materials that exhibit non-zero anisotropy. These results demonstrate that, when\nprovided with an appropriately anisotropy-rich dataset, our model can reliably learn the\nunderlying symmetry relationships encoded in the crystal structure and map them to the\ncorresponding tensorial spectra. The auxiliary dataset and the associated model checkpoints have\nalso been made publicly available in our GitHub repository.\nResponse Figure 1 (Supplementary Figure 16). Distribution of full-tensor MAE values across the test\nset of the auxiliary dataset (left) and representative spectra obtained by stratified sampling across quartiles\n(right). From the cumulative kernel-density-estimator (KDE) plot, we randomly selected two systems in\neach quartile (𝑄 −𝑄 ), resulting in 8 representative examples that span the full error distribution.\n! \"\nResponse Figure 2 (Supplementary Figure 17). Distribution of full-tensor MAE values across the test\nset of the auxiliary dataset (left) and representative spectra obtained by stratified sampling across quartiles\n(right). From the cumulative kernel-density-estimator (KDE) plot, we randomly selected two systems in\neach quartile (𝑄 −𝑄 ), resulting in 8 representative examples that span the full error distribution.\n! \"\nSummary of revisions made based on preceding discussion:\n• Main text, page 13, line 352—358.\n• Supplementary Section 8. Line 411 – 438.\n• Supplementary Table 6.\n• Supplementary Figure 16.\n• Supplementary Figure 17.\nThe whole paper and Supplemental Information are quite technical with a lot of definitions\nrepeated for convenience of the reader from original sources.\nWe thank the reviewer for the observation. Given the technical nature of the work, we\nintentionally included brief definitions from standard sources to keep the exposition self-\ncontained and accessible to a broad readership. We hope this improves clarity for the broad\nreadership of the journal readers who may not be familiar with some aspects of the underlying\nformalism.\nIf am not mistaken, the statement \"The dielectric tensor can be expressed as the outer product of\ntwo rank-1 vectors\" on page 3, rephrased as Eq. (1), is simply incorrect. Any tensor can trivially\nbe written as a sum of outer products.\nFurthermore, discussing these outer products is superfluous because it only leads to equations (3)\nthrough (5), which describe the action of the e3nn package, i.e., Ref. [21]. I strongly recommend\nstreamlining or fixing this part before publishing the manuscript.\nIf I misunderstood something, the authors should try to prevent that others make the same\nmistake.\nWe thank the reviewer for pointing out the ambiguity in our original phrasing. We agree that a\ngeneral dielectric tensor cannot be expressed as a single outer product of two vectors, and that\nour wording could be misinterpreted as implying a rank-1 structure. In response, we have\nremoved this part of the main text, as it was superfluous, and now refer the reader to\nSupplementary Section 6 for a detailed derivation. We have revised the text accordingly to\neliminate this ambiguity and to streamline the connection to the spherical harmonic\nrepresentation, thereby shortening this part of the discussion.\nSummary of revisions made based on preceding discussion:\n• Main text, page 3, line 105— page 4, line 110.\nA very minor suggestions is to define \\sigma_1,2 in Eq. (7) with a sentence like “where \\sigma_1\nand \\sigma 2 are learnable parameters….” [Taken from page 7 of the SI].\nWe thank the reviewer for the helpful suggestion. We have added the corresponding changes for\n𝜎 and 𝜎 .\n* +\nSummary of revisions made based on preceding discussion:\n• Main text, page 6, line 173.\nIncidentally, the respective equation on page 7 of the SI is unnumbered, unlike almost all other\nequations. On purpose?\nRegarding the Supplemental Information, we thank the reviewer for noting the inconsistency in\nequation numbering; this was unintentional, and we have corrected it in the revised SI.\nSummary of revisions made based on preceding discussion:\n• Supplementary Information, page 7, line 142.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nWe sincerely thank Reviewer #2 for their thoughtful evaluation and for participating in the\nNature Communications co-review initiative. We appreciate the time and care invested in\nassessing our work.\nReviewer #3 (Remarks to the Author):\nThe authors have addressed my concerns, and thus I support the publication.\nWe thank Reviewer #3 for their positive feedback and for acknowledging that our revisions have\naddressed the earlier concerns. We appreciate the support for publication.\nReviewer #4 (Remarks to the Author):\nWe thank the authors taking their time to revise the manuscript according to reviewers'\ncomments. We believe that the manuscript deserves publication in its current form.\nWe sincerely thank Reviewer #4 for their careful reading and for recognizing the improvements\nmade in our revision. We appreciate the supportive assessment of the manuscript in its current\nform.\nReviewer #5 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nWe thank Reviewer #5 for their contribution through the Nature Communications co-review\ninitiative and for taking the time to evaluate our revised manuscript. We greatly appreciate the\nfeedback provided.\nSincerely,\nQimin Yan\non behalf of all co-authors\nDepartment of Physics and\nQuantum Materials and Sensing Institute\nNortheastern University"}, {"round": 3, "reviewer_comments": {"reviewer_1": "I am pleased with the authors' constructive response to our critical comments. The new data and figures on low-symmetry\ncompounds shown in the SI will certainly help readers to assess strengths and weaknesses of the work and to expand the\nmodel if necessary.\nI can therefore recommend the manuscript for publication without further changes.\n(Remarks on code availability)", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.\n(Remarks on code availability)"}, "author_response": ""}]

Decision letter