Mechanism of glass transition temperature enhancement in multicomponent epoxy resins incorporating triazine rings

Abstract

Epoxy resins are widely used in electronic and structural applications. However, their use at high temperatures is limited by their glass transition temperature (T_g). Although the addition of multifunctional resins increases T_g via crosslink density, the molecular-scale contribution of specific intermolecular interactions, such as π–π stacking, remains poorly understood in multicomponent systems. Here, we investigated the molecular-level mechanism of T_g enhancement in a multicomponent epoxy resin system incorporating a triazine-based epoxy resin, tris(2-epoxypropyl)isocyanurate (TEPIC), into a conventional epoxy resin. The experimental results revealed a non-monotonic dependence of T_g on TEPIC content, with an initial decrease at low concentration followed by a pronounced increase at higher concentration, which cannot be explained solely by changes in crosslink density. Structural analysis using wide-angle X-ray scattering and molecular dynamics (MD) simulations showed that the amorphous halo originating from intermolecular packing splits into two peaks, suggesting that one of these peaks includes a contribution from π–π stacking interactions. MD simulations showed that ring pairs with centre-of-mass distances of 3.4–5.8 Å form stable π–π stacking structures, and that the population of these interactions correlates strongly with the observed variation in T_g. Quantum chemical calculations further demonstrated that benzene–triazine stacking interactions introduced by TEPIC are significantly stronger than benzene–benzene stacking. These results indicate that the enhancement of T_g arises not only from an increase in crosslink density but also from localized intermolecular constraints induced by strong stacking interactions, providing a molecular-level design guideline for high-performance epoxy resins.

---

1. Introduction

Epoxy resins are widely used in advanced electronic materials, such as adhesives and semiconductor encapsulants, owing to their high mechanical strength and excellent heat resistance, electrical insulation, and adhesion properties. In recent years, the continual miniaturization of electronic devices and the increasing demand for high-performance devices have led to a significant increase in power density and heat generation inside these devices. Therefore, the resin materials used for insulation, protection, and mechanical support in electronic devices must operate under high-temperature conditions. Epoxy resins exhibit a marked degradation in their thermal and mechanical properties at temperatures above their glass transition temperature (T_g),^1–3 so they are usually limited to application temperatures below T_g. Therefore, it is necessary to increase the T_g of epoxy resins to enable further improvements in the performance of miniaturized devices.

Diglycidyl ether of bisphenol A (DGEBA) is widely used owing to its relatively low cost and ease of handling; however, its relatively low T_g limits its use under demanding thermal conditions. To overcome this drawback, DGEBA is often modified by incorporating multifunctional epoxy resins to increase its crosslink density and hence, the T_g of the composite system.^4–9 Despite the wide range of strategies proposed for increasing T_g, the fundamental details of the resulting crosslinked network structures remain insufficiently understood because the direct experimental observation of epoxy networks at the molecular level is inherently challenging. Therefore, many studies have indirectly inferred changes in the internal structure from macroscopic properties, such as the T_g and elastic modulus.

X-ray scattering methods are effective techniques for evaluating the periodicity of crosslinked structures and intermolecular interactions within epoxy resins.^10–15 Ogawa et al.¹⁰ investigated the structural periodicity of epoxy resins using medium-angle X-ray scattering and revealed that a periodic structure at q = 1.2 nm⁻¹ centered around the curing agent is formed as the epoxy–amine curing reaction proceeds. Dwyer et al.¹¹ analyzed the polymer network structure of epoxy resins by wide-angle X-ray scattering (WAXS) and reported that the distance between π–π stacking rings increases upon heating.

To obtain a more detailed understanding of the relationship between the crosslinked structure and macroscopic properties, structural analysis at the molecular level is essential. Molecular dynamics (MD) simulations are a powerful technique for reproducing the curing process and network structure of epoxy resins at the molecular scale. In recent years, MD techniques have been increasingly used in combination with experimental techniques such as WAXS.^10,16–22 Li et al.²³ investigated the relationship between π–π stacking and the T_g of uncured silicon-containing arylacetylene resins using MD simulations. They demonstrated that increasing the fraction of aromatic rings enhances π–π interactions, which in turn contribute to increasing the T_g and viscosity. Yamamoto et al.¹⁵ investigated the correlation between crosslinked structures and physical properties in epoxy resins with various stoichiometric ratios using WAXS and MD simulations. They reported that, as the amine content increases beyond the stoichiometric ratio, the expansion of the distance between phenyl groups during the curing reaction is suppressed, providing steric resistance during deformation, leading to an increase in Young's modulus. Kawagoe et al.¹⁷ calculated scattering intensities using MD simulations and showed, through comparison with WAXS measurements, that the amorphous halos and peaks emerging during curing were successfully reproduced by the simulations. Furthermore, Zhao et al.⁸ used MD simulations to demonstrate that a higher tetraglycidyl diamino diphenyl methane (TGDDM) content in the DGEBA/TGDDM/4,4′-diaminodiphenyl sulfone (4,4′-DDS) system increases the number of ring structures, thereby improving the thermomechanical properties.

Most previous studies primarily interpreted the thermomechanical properties of cured epoxy resins in terms of network formation, via parameters such as the crosslink density and segmental rigidity (e.g., aromatic versus aliphatic linkages between junctions). In contrast, fewer studies have leveraged the richer structural information accessible by WAXS, including short-range intermolecular ordering, and combined it with MD simulations to establish structure–property relationships beyond network connectivity. Therefore, to deepen the understanding of structure–property relationships beyond network connectivity, we investigated a multicomponent epoxy resin system in which a triazine-based epoxy resin, tris(2-epoxypropyl)isocyanurate (TEPIC), is incorporated into a conventional DGEBA matrix. TEPIC has been widely used as an additive to enhance thermomechanical properties.^24–26 In this study, the internal structure and thermomechanical properties were predicted by MD simulations, and the validity of the simulations was examined through comparison with WAXS measurements. This study aimed to elucidate the mechanism of T_g enhancement induced by the introduction of triazine rings based on the correlation between the internal structure and thermomechanical properties.

2. Experiments

2.1. Preparation

DGEBA (Mitsubishi Chemical Corporation) and TEPIC (Nissan Chemical Corporation) epoxy resins and 4,4′-DDS (Tokyo Chemical Industry Co., Ltd) curing agent were used. The molecular structures of these materials are shown in Fig. 1. Specimens were prepared by varying the TEPIC content in the DGEBA/TEPIC/4,4′-DDS resin system (0, 5, 10, 20, and 30 wt%), which were experimentally analyzed using thermomechanical and WAXS techniques.

Fig. 1 Structures of (a) base resins and (b) curing agent.

The cured epoxy resins were prepared as follows. First, DGEBA and 4,4′-DDS were stirred at 120 °C for 45 min. Subsequently, TEPIC preheated to 130 °C was added, and the mixture was stirred for an additional 10 min under vacuum using a rotation–revolution vacuum mixer (ARV-310P, Thinky) to obtain a homogeneous resin mixture. The 0 wt% mixture was prepared using the exact same procedure, simply omitting TEPIC. For all component formulations, the mixtures were carefully formulated to maintain a stoichiometric ratio of 1.0 between the total number of epoxy groups (from DGEBA and TEPIC) and the amine groups (from 4,4′-DDS). The mixture was poured into a mold made of glass plates covered with release film, and cured in an oven by heating from room temperature to 180 °C at a rate of 1.7 °C min⁻¹, followed by holding at 180 °C for 2 h. The specimen thickness was controlled to 2 mm using polytetrafluoroethylene spacers. Various test specimens were prepared by cutting the cured plates to the required dimensions.

2.2. Thermomechanical properties

Two methods, thermomechanical analysis (TMA) (TMA7300, Hitachi High-Tech) and dynamic mechanical analysis (DMA) (DMA7100, Hitachi High-Tech), were used to determine the T_g values of the cured resins. In the TMA measurements, the specimens were heated from room temperature to a temperature above the T_g at a rate of 2 °C min⁻¹ under compression mode with an applied load of 1 mN. The specimen size was 10 mm × 10 mm × 2 mm, and T_g was determined from the intersection of the linear fits to the temperature–thermal strain relationship in the glassy and rubbery regions along the 10 mm length direction. For DMA, the specimens were heated from room temperature to a temperature above T_g, and T_g was evaluated from the peak position of the loss modulus (E″). E″ was selected over the storage modulus (E′) and loss tangent because it provided superior sensitivity and resolution for distinguishing the multiple glass transitions observed in the multicomponent systems. The measurements were conducted in tensile mode with a heating rate of 2 °C min⁻¹ and a frequency of 1 Hz. The specimen dimensions were 30 mm × 10 mm × 2 mm. The gauge length was set to 15 mm, and the displacement was set to 7.5 µm.

2.3. Wide-angle X-ray scattering measurements

WAXS measurements were employed to quantitatively evaluate the molecular structure inside the resin. WAXS measurements were conducted using BL08W at the NanoTerasu synchrotron facility (Miyagi, Japan). 2D-WAXS patterns were recorded using a 2D hybrid photon counting detector Eiger2 R (Bruker K.K., Japan). WAXS measurements were performed using a camera length of 80 mm at room temperature and over the range of 80–280 °C in 20 °C increments.

3. Numerical methods

3.1. Quantum-chemical reaction-path calculations

In the DGEBA/TEPIC/4,4′–DDS resin system, the epoxy–amine network can be formed through six elementary reaction types: (i) DGEBA reacting with primary amines; (ii) DGEBA reacting with secondary amines generated after DGEBA addition (DGEBA/secondary amine/DGEBA); (iii) DGEBA reacting with secondary amines generated after TEPIC addition (DGEBA/secondary amine/TEPIC); (iv) TEPIC reacting with primary amines; (v) TEPIC reacting with secondary amines generated after DGEBA addition (TEPIC/secondary amine/DGEBA); and (vi) TEPIC reacting with secondary amines generated after TEPIC addition (TEPIC/secondary amine/TEPIC). However, because these six reaction types occur simultaneously during curing, it is difficult to experimentally evaluate the relative magnitudes of the activation energies associated with each reaction. Therefore, only the reactions that could be evaluated by differential scanning calorimetry (DSC) measurements were considered in this study, i.e., DGEBA/primary amine, DGEBA/secondary amine/DGEBA, TEPIC/primary amine, and TEPIC/secondary amine/TEPIC. The activation energies and heats of formation, which characterize the interactions between the base resin and curing agent, are key factors governing both the curing process and the properties of the cured resin.^27,28 These reaction parameters were determined using the global reaction route mapping (GRRM) method^29–32 combined with density functional theory (DFT) calculations at the B3LYP/6-31G(d) level.

Table 1 summarizes the activation energies and heats of formation obtained from the GRRM calculations. Based on the order of the activation energies, the reaction between TEPIC and primary amines proceeds first, followed by reactions between DGEBA and primary amines and between DGEBA and secondary amines, and finally the reaction between TEPIC and secondary amines occurs. This reaction sequence is in good agreement with the results obtained from separately conducted DSC measurements. These parameters were used as input data for the MD simulations of epoxy curing. For the DGEBA/secondary amine/TEPIC reaction, the reaction heat and activation energy of the DGEBA/secondary amine/DGEBA reaction were used, while for the TEPIC/secondary amine/DGEBA reaction, the corresponding values of the TEPIC/secondary amine/TEPIC reaction were adopted. The values for the DGEBA/primary amine and DGEBA/secondary amine/DGEBA reactions were taken from previous studies.⁸

Table 1 Activation energies E_a and heat of formation H_f obtained by GRRM

Reaction type	Ea [kcal mol^−1]	Hf [kcal mol^−1]
DGEBA/primary amine^8	50.01	25.04
DGEBA/secondary amine/DGEBA^8	41.10	20.33
TEPIC/primary amine	48.58	18.35
TEPIC/secondary amine/TEPIC	54.90	11.35

---

3.2. Molecular dynamics simulations

3.2.1. Curing simulation. A curing MD simulation framework developed by Okabe et al.³³ effectively predicted the crosslinked structures and properties of epoxy resins formed through chemical reactions. This framework employs a reaction model derived from quantum chemical calculations and has been successfully applied to a wide range of resin systems, demonstrating high accuracy in predicting thermomechanical properties.^8,17,34–39 In this study, DGEBA, TEPIC, and 4,4′-DDS molecules were randomly placed in a simulation box of 200 Å × 200 Å × 200 Å at 300 K, with the total number of atoms set to approximately 40 000. Periodic boundary conditions were imposed in all directions, and the simulation box was isothermally compressed at 300 K until the system density reached 0.9 g cm⁻³. Subsequently, a 300 ps NPT simulation was performed at 300 K and 1 atm to obtain the initial uncured structure.

Curing MD simulations were then performed by applying the following reaction model to the initial structure. Two criteria were applied to the reaction model. The first criterion was a distance condition, which determines candidate reactive pairs. When the distance between the representative atoms of two functional groups falls below the reaction cutoff distance R_c, the pair is considered a reaction candidate. In this work, the terminal C atom of the epoxy group and the N atom of the amine group were defined as the representative atoms, and the cutoff distance was set to R_c = 5.64 Å. The second criterion was the reaction probability condition, which determines whether a candidate pair will react. The reaction probability p for pairs satisfying the distance condition was calculated using the Arrhenius equation:

here, E_a is the activation energy obtained from the GRRM calculations, R is the universal gas constant, T is the local temperature at the reactive site obtained from the MD simulation, and A is an acceleration factor, which was set to 10¹⁴ (ref. 40) to allow the curing reactions to be completed within a realistic computational time.

If the calculated p exceeds a uniformly generated random number a in the range from 0 to 1, a crosslinking reaction is assumed to occur between the corresponding functional groups. After each reaction event, structural relaxation was performed using the conjugate gradient method, followed by an NPT simulation (15 ps, 1 atm) and an NVT simulation (1 ps). Finally, the atomic velocities in the reaction region were scaled so that the kinetic energy satisfied K_after = K_before + H_f, where K_before and K_after are the total kinetic energies of the reactive atoms before and after the reaction, respectively, and H_f is the heat of formation obtained from the GRRM. The curing simulation was repeated until no further reactions occurred.

Before evaluating the thermomechanical properties of the cured models obtained from the curing simulations, annealing was performed following the methodology reported in previous studies.^35,40 The annealing process consisted of an NVT simulation in which the temperature was cyclically varied between 300 and 1000 K, followed by an NPT simulation for equilibration. All MD simulations in this study were performed using the LAMMPS MD simulator,⁴¹ and the crosslink-formation procedure was implemented using in-house code. Interatomic interactions were described using the DREIDING force field,⁴² and partial charges were assigned using the charge equilibrium (QEq) method⁴³ at each cycle. The QEq parameters reported by Zhang et al.⁴⁴ and Demir and Walsh⁴⁵ were employed. The van der Waals interactions were modeled using a Lennard–Jones potential with a cutoff distance of 12 Å, and electrostatic interactions were calculated using the particle–particle–particle–mesh method.⁴⁶ The r-RESPA algorithm⁴⁷ was used for time integration, with different timesteps assigned to specific interactions: 0.2 fs for bond and angle interactions, and 1.0 fs for dihedral, van der Waals, and electrostatic interactions. For both NVT and NPT simulations, temperature and pressure were controlled using the Nosé–Hoover thermostat and barostat, respectively.⁴⁸

3.2.2. Evaluation of glass transition temperature. The cured resin models after annealing were used to evaluate T_g, which was determined from stepwise cooling simulations in the NPT ensemble. The simulation box was first heated to 700 K and subsequently cooled in steps of 20 K to 200 K. At each temperature step, the system was equilibrated for 500 ps, after which the volume was recorded. The T_g was defined as the intersection of the linear fits to the temperature-volume relationship in the glassy and rubbery regions.

3.2.3. Evaluation of structural properties. To evaluate the atomic-scale structure, the partial radial distribution function g_αβ(r) was calculated for the annealed structure (300 K) ashere, N_αβ(r) denotes the number of α–β atomic pairs within the distance range from r to r + dr, and n_β is the number density of species β.

The coherent scattering intensity was then evaluated from the Fourier transformation of g_αβ(r). Because the transformation was performed over a finite range (0 ≤ r ≤ r_max), a window function was applied to suppress termination ripples. In this study, the Lorch window function,

was employed, where r_max was set to half of the MD simulation cell length. Using this window function, the coherent scattering intensity I_coh was calculated from the partial radial distribution functions ashere, N is the total number of atoms, c_i = N_i/N is the atomic fraction of the ith species, f_i(q) is the atomic scattering factor of the ith species in X-ray scattering, and n₀ denotes the average total number density. In addition, q= 4π sin θ/λ is the magnitude of the scattering vector, where λ is the wavelength of the incident X-ray radiation and 2θ is the diffraction angle.

3.3. Quantum-chemical evaluation of fragment interaction energies

Quantum-chemical calculations were performed using the Gaussian program package to evaluate interaction energies between selected fragment pairs. Geometry optimization steps were performed at the MP2/6-311G(d,p) level. Single-point interaction energies were then computed at the MP2/aug-cc-pVDZ level with basis set superposition error correction using the counterpoise method,⁴⁹ and the resulting interaction energies are presented in the subsequent sections.

4. Results and discussion

4.1. Thermal properties

Fig. 2 shows the results of the DMA and TMA measurements and summarizes the T_g values for each resin determined from these measurements and MD simulations. The DMA results show that, compared with the unmodified system (0 wt%), adding 5 wt% TEPIC decreases T_g, which then increases markedly when the TEPIC content exceeds 10 wt%. Similarly, the TMA results indicate lower T_g values for the 5 wt% and 10 wt% TEPIC systems compared to that of the unmodified system, whereas a pronounced increase in T_g is observed upon the addition of 20 wt% or more TEPIC. Thus, both DMA and TMA measurements consistently reveal a characteristic behavior in which a decrease in T_g is observed with the addition of low TEPIC concentrations. This trend is also supported by the T_g values obtained from MD simulations. In contrast, an examination of E′ measured by DMA at temperatures above the T_g shows that the crosslink density increases monotonically with increasing TEPIC content. In general, the incorporation of resins with higher density of functional groups increases the crosslink density, thereby increasing T_g;²⁶ however, this effect alone cannot sufficiently explain the observed decrease in T_g at low TEPIC additions. Supporting this interpretation, Downey et al.⁵⁰ showed that adding an aliphatic trifunctional epoxy without benzene rings to a DGEBA/amine curing agent system led to a decrease in T_g with increasing additive content despite increased crosslinking. Therefore, the mechanism underlying this behavior was further investigated by analyzing temperature-dependent WAXS patterns and their correlation with the simulated internal structures.

Fig. 2 Thermal properties of each resin. (a) Temperature dependence of the storage modulus (E′) and loss modulus (E″) measured by DMA. Markers on the E″ indicate T_g. (b) Temperature dependence of the thermal strain obtained by TMA. Markers indicate T_g. (c) T_g as a function of TEPIC concentration measured by DMA, TMA, and MD simulations. For the 0 wt% sample, TMA data measured by Odagiri et al.³⁶ are shown.

Furthermore, the DMA results indicate that some compositions exhibit two T_g values. In particular, the system containing 20 wt% TEPIC shows the coexistence of clearly distinguishable rigid and soft regions. Importantly, this behaviour is unlikely to originate from mesoscopic phase separation into DGEBA–rich and TEPIC–rich phases: despite the TEPIC content being at most 30 wt%, the lower-temperature peak significantly weakens when the TEPIC content is increased from 20 to 30 wt%, rather than persisting as would be expected for stable phase coexistence. Instead, the results are consistent with network heterogeneity, i.e., coexistence of (i) three-component regions (DGEBA/TEPIC/4,4′-DDS) where TEPIC effectively contributes to local motion restriction and (ii) two-component regions (DGEBA/4,4′-DDS) where the constraint effect is weaker. This interpretation is further supported by the MD results in the following section, which show the same qualitative tendency without indications of large-scale demixing. Furthermore, separate small-angle X-ray scattering measurements showed no scattering peaks indicative of phase separation.

4.2. Internal structure indicated by WAXS

Fig. 3 shows experimental and MD-simulated WAXS patterns. For all resin systems, the experimental WAXS patterns show four major peaks, which are henceforth called peak 1, peak 2, peak 3, and peak 4, in order from the low-angle side. Fig. 3(c) shows that the MD simulations successfully reproduced three of the four peaks (peaks 2–4) observed in the experimental WAXS patterns. In addition, the simulations captured both the temperature-induced shift of peak 2 and the attenuation of the shoulder with increasing temperature. The good agreement between the experimental and simulated patterns indicates the reliability of the simulated structures for interpreting the molecular origin of peaks 2–4. The absence of peak 1 in the simulated patterns is likely due to the limited simulation box size (approximately 70 × 70 × 70 Å), which is insufficient to represent the larger-scale correlations (approximately 18 Å) that give rise to low-angle scattering. Previous studies have reported that peak 4 originates from the average C–C distance,¹¹ while peaks 2 and 3 correspond to the amorphous halo originating from intermolecular correlations^13,17,51 or π–π stacking.¹¹ For the material system employed in this study, peak 1 is thought to originate from the distance between the sulfur atoms of 4,4′-DDS molecules connected by a DGEBA unit (see SI Section S1).

Fig. 3 X-ray scattering intensities obtained from WAXS measurements and MD simulations. (a) WAXS profiles at 33 °C. (b) Temperature-dependent WAXS profiles. (c) Scattering profiles calculated from MD simulations for 0 wt% and 30 wt% TEPIC, which are vertically shifted to avoid overlap.

Peaks 2 and 3 overlap, making it difficult to directly observe how each peak responds to changes in composition or temperature. Therefore, in this study, peak separation was performed using Gaussian functions in the amorphous halo region, allowing the two peaks to be deconvoluted, as shown in Fig. 4. Fig. 5 shows the temperature dependence of the structural periods calculated from the position corresponding to the peak maximum as d = 2π/q for peaks 2 and 3. Both peaks 2 and 3 exhibit an increase in the corresponding structural periodicity with increasing temperature, and the slope of the temperature dependence shows an inflection point near T_g. This is in good agreement with the linear thermal expansion behavior observed from the TMA measurements, suggesting that these peaks reflect structural information originating from intermolecular correlations. Fig. 6 shows the temperature dependence of the area ratio of peak 3 relative to the total area of the amorphous halo. The area of peak 3 normalized by the total amorphous halo area decreases with increasing temperature. This decrease becomes more pronounced around T_g. This suggests that some of the structures contributing to peak 3 may weaken as molecular mobility increases at higher temperatures.

Fig. 4 Peak separation of the scattering profile for the 0 wt% TEPIC sample at (a) 33 °C and (b) 280 °C. The profiles obtained from experimental WAXS measurements were fitted using Gaussian functions.

Fig. 5 Temperature dependence of the structural periodicity of (a) peak 2 and (b) peak 3 derived from experimental WAXS profiles. The structural period values were calculated from the position of each peak top.

Fig. 6 Temperature dependence of the area ratio of peak 3 relative to the total area of the amorphous halo derived from experimental WAXS profiles. The amorphous halo area is defined as the sum of the integrated areas of peaks 2 and 3 obtained by Gaussian fitting.

4.3. π–π stacking

The WAXS results suggest the presence of ordered structures that weaken with increasing temperature. In epoxy systems containing aromatic rings, such thermally unstable ordered structures are often associated with π–π interactions. Therefore, MD simulations were used to examine the presence and quantify the number of π–π stacking pairs for each composition. The structural reliability of the MD models was verified in the previous section through comparison with experimental WAXS. In the present resin system, π–π interactions can occur between benzene rings and between electron-rich benzene rings and electron-deficient triazine rings.^52–54

For each composition, NPT simulations were performed for 400 ps and trajectories were recorded every 10 ps. From the obtained trajectories, ring pairs were considered to be in a π–π stacked state if they satisfy the criteria of a center-of-mass distance d ≤ 6 Å and an angle α ≤ 30° between the normal vectors of the two ring planes (see Fig. 7). Only pairs that met these criteria throughout the entire 400 ps trajectory were counted as stacked pairs. The analysis results were averaged over N = 10 independent simulations.

Fig. 7 Number of π–π stacked rings as a function of TEPIC content, classified into total, benzene–benzene, and benzene–triazine. The number of rings was normalized by the total number of rings within the simulation box to enable comparisons between systems with different compositional ratios. The inset shows a schematic of two aromatic rings and the geometric criteria of their relative position and orientation used to define π–π stacking.

Fig. 7 shows the total number of stacking pairs, as well as the numbers of benzene–benzene and benzene–triazine stacking pairs, for each TEPIC content. The number of pairs was normalized by the total number of rings within the simulation box to enable the comparison of systems with different compositional ratios. Fig. 8 shows simulation boxes in which only π–π stacked ring structures are shown for each TEPIC content.

Fig. 8 MD simulation boxes showing only π–π stacked rings. Gray: C atoms; blue: N atoms.

For compositions containing ≥10 wt% TEPIC, both benzene–triazine and benzene–benzene stacking were observed. In contrast, triazine–triazine stacking was not identified based on our criteria for any of the compositions investigated. Interestingly, although the concentration of benzene rings in the system decreases with increasing TEPIC content, the number of benzene–benzene stacking pairs does not decrease. This behavior is attributed to the addition of the relatively short-chain TEPIC molecules, which facilitate configurations in which benzene rings can approach each other more closely, thereby promoting stacking interactions. Fig. 9 shows the distribution of the center-of-mass distances for stacked ring pairs. The stacked-pair center-of-mass distances span a range of 3.4–5.8 Å, which is comparable to the real-space correlation lengths derived from peak 3.

Fig. 9 Histograms of center-of-mass distances for the identified π–π stacking pairs: for (a) 0 wt% and (b) 30 wt% TEPIC systems.

Furthermore, the composition dependence of the total number of stacking pairs reveals an initial decrease from 0 wt% to 5 wt%, followed by a monotonic increase for TEPIC contents of ≥10 wt%. This trend agrees well with the composition dependence of T_g, suggesting that, in addition to network connectivity (crosslink density), π–π stacking provides local intermolecular constraints that contribute to the observed behavior of the T_g of the epoxy network. The minimum T_g observed at 5 wt% TEPIC addition is attributed to the low TEPIC content, which is insufficient to form benzene–triazine stacking or promote closer approaches between benzene rings. Therefore, the introduction of 5 wt% TEPIC simply reduces the number of benzene–benzene stacking interactions.

4.4. Evaluation of π–π interaction energies

Finally, quantum-chemical calculations were performed for benzene–benzene and benzene–triazine ring fragments to evaluate their π–π stacking interaction energies. In Fig. S3 of the SI, benzene–benzene stacking adopts a slipped configuration, whereas benzene–triazine stacking forms a face-to-face configuration. Table 2 summarizes the interaction energies, center-of-mass distances, and configuration for each stacking structure. Benzene–triazine stacking exhibits a stronger interaction energy than benzene–benzene stacking, suggesting that benzene–triazine stacking is more resistant to thermal disruption and more effectively constrains the molecular motion of rings. The stacking distances between ring pairs obtained from quantum chemical calculations are shorter than those obtained from MD simulations. This is attributed to steric hindrance induced by the crosslinked structures in the MD simulations, which prevents the rings from approaching as closely as they would intrinsically prefer.¹⁵

Table 2 Calculated stacking energies, center-of-mass (CoM) distances, and preferred configurations for different ring pairs

Type	Energy [kcal mol^−1]	CoM distance [Å]	Configuration
Benzene–benzene	−4.03	3.38	Slip
Benzene–triazine	−7.26	3.28	Face-to-face

---

The MD simulation results indicate that a higher number of π–π stacking pairs correlates with an increase in T_g. Furthermore, the quantum chemical calculations reveal that benzene–triazine stacking introduced by TEPIC addition suppresses molecular motion more strongly than benzene–benzene stacking. Taken together, these findings suggest that the increase in T_g upon sufficient TEPIC addition is not governed solely by crosslink density, but also by strong, localized intermolecular interactions associated with ring stacking within the network structure.

5. Conclusions

A multicomponent epoxy resin system was investigated in which a triazine-based epoxy resin, TEPIC, was added to a DGEBA/4,4′-DDS system. By combining WAXS, thermomechanical analysis, and MD simulations, the relationship between the internal structure and the glass transition behavior was systematically examined.

Experimental thermal characterization revealed a characteristic decrease in T_g with the addition of a low TEPIC concentration, followed by a steady increase in T_g at higher TEPIC contents. This non-monotonic behavior cannot be explained solely by an increase in crosslink density induced by the addition of TEPIC.

The origin of the peaks observed in the WAXS measurements was clarified by comparison with MD simulation results. The amorphous halo originating from intermolecular packing was found to split into two peaks, suggesting that one of these peaks includes a contribution from π–π stacking interactions. Detailed MD simulations revealed that a subset of ring pairs with center-of-mass distances in the range of 3.4–5.8 Å participates in π–π stacking. Furthermore, the total number of such stacking pairs showed a trend consistent with the observed variation in T_g. Accordingly, π–π stacking interactions play a key role in determining the glass transition behavior of the epoxy network.

At high TEPIC contents, stacking between benzene and triazine rings occurs, and quantum chemical calculations revealed that benzene–triazine stacking exhibits stronger stacking energies than benzene–benzene stacking. These strong stacking interactions remain relatively stable against thermal fluctuations and strongly restrict the motion of rings. Such suppression of molecular mobility throughout the network induced by stacking interactions is considered to be responsible for the increase in T_g observed at high TEPIC contents.

In conclusion, the mechanism underlying the enhancement of T_g upon the introduction of triazine rings is attributed to both an increase in crosslink density and the introduction of localized intermolecular constraints arising from stacking interactions. This study presents a molecular-scale mechanism that rationally explains the experimentally observed characteristic T_g behavior and clarifies the role of π–π interactions in determining the thermal properties of epoxy resins. These findings provide a molecular-level design guideline for high-performance epoxy resins.

Author contributions

Yuuki Kinugawa: writing – original draft, writing – review & editing, validation, investigation, formal analysis. Yoshiaki Kawagoe: writing – review & editing, methodology, validation, investigation, formal analysis. Keigo Matsumoto: writing – review & editing, validation, investigation. Masashi Ohno: writing – review & editing, validation, investigation. Shoko Mishima: writing – review & editing, validation, investigation. Takahiko Kawai: writing – review & editing, validation, investigation. Tomonaga Okabe: writing – review & editing, methodology, funding acquisition, supervision.

Conflicts of interest

There are no conflicts to declare.

Data availability

The code for LAMMPS and curing python code can be found at https://www.lammps.org/#gsc.tab=0 and http://www.cosmic.plum.mech.tohoku.ac.jp/english.html, respectively.

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6cp01039c.

Acknowledgements

This study was partially supported by JST, K Program (Grant Number JPMJKP24W1), and by JSPS KAKENHI (Grant Number 23K25987). Numerical simulations were performed using the AFI-NITY supercomputer system at the Advanced Fluid Information Research Center, Institute of Fluid Science, and the AOBA supercomputer system at the Cyberscience Center, Tohoku University, Japan. We are grateful to the staff of QST NanoTerasu Center and PhoSIC for their support during our experiments.

References

1. R. Rohini and S. Bose, Phys. Chem. Chem. Phys., 2015, 17, 7907–7913.
2. A. T. Sunny, R. Adhikari, S. Mathew and S. Thomas, et al., Phys. Chem. Chem. Phys., 2016, 18, 19655–19667.
3. K. Ryuzono, S. Koo, Y. Hoshikawa, Y. Kawagoe and T. Okabe, Composites, Part A, 2025, 196, 108963.
4. A. Bajpai, J. R. Davidson and C. Robert, Appl. Mech., 2021, 2, 419–430.
5. G. Cicala, S. Mannino, A. Latteri, G. Ognibene and G. Saccullo, Adv. Polym. Technol., 2018, 37, 1868–1877.
6. J. Anurangi, S. Jayalath, S. R. Senevirathna, M. Herath and J. Epaarachchi, Polym. Compos., 2024, 45, 10898–10910.
7. D. J. Hourston, J. M. Lane and H. X. Zhang, Polym. Int., 1997, 42, 349–355.
8. Y. Zhao, G. Kikugawa, K. Shirasu, Y. Kawagoe and T. Okabe, Polymer, 2024, 297, 126817.
9. X. Zhang, Q. Chi, C. Tang, H. Li, C. Zhang, Z. Li and T. Zhang, J. Mater. Sci.: Mater. Electron., 2023, 34, 638.
10. H. Ogawa, M. Aoki, S. Ono, Y. Watanabe, S. Yamamoto, K. Tanaka and M. Takenaka, Langmuir, 2022, 38, 11432–11439.
11. D. B. Dwyer, S. Isbill, Z. E. Brubaker, J. K. Keum, W. Bras and J. L. Niedziela, ACS Appl. Polym. Mater., 2023, 5, 5961–5971.
12. J. Liu, X. Liu, X. Cui, J. Qin, M. Wu, L. Fu, M. Shi, D. Wang and L. Liang, Polym. Chem., 2023, 14, 1665–1679.
13. S. Kumar and W. W. Adams, Polymer, 1987, 28, 1497–1504.
14. R. Lovell and A. Windle, Polymer, 1990, 31, 593–601.
15. S. Yamamoto, N. T. Phan, K. Kihara, A. Shundo and K. Tanaka, Polym. J., 2025, 57, 357–366.
16. Q. Yang, X. Yang, X. Li, L. Shi and G. Sui, RSC Adv., 2013, 3, 7452–7459.
17. Y. Kawagoe, Y. Kinugawa, K. Matsumoto, M. Ohno, N. Kishimoto, T. Kawai and T. Okabe, Phys. Chem. Chem. Phys., 2024, 26, 24250–24260.
18. K. Hackenstrass, N. T. Jonasson, M. Hartwig-Nair, T. Rosén, S. Florisson and M. Wohlert, Faraday Discuss., 2026, 263, 52–64.
19. R. Huo, Z. Zhang, N. Athir, Y. Fan, J. Liu and L. Shi, Phys. Chem. Chem. Phys., 2020, 22, 19735–19745.
20. K. Yamaguchi, D. Kawaguchi, N. Miyata, T. Miyazaki, H. Aoki, S. Yamamoto and K. Tanaka, Phys. Chem. Chem. Phys., 2022, 24, 21578–21582.
21. J. Kang, C. Wang, D. Li, G. He and H. Tan, Phys. Chem. Chem. Phys., 2015, 17, 16519–16524.
22. J. Zhu, L. Wang, J. Lin, L. Du and Q. Zhuang, Phys. Chem. Chem. Phys., 2021, 23, 14027–14036.
23. H. Li and W. Zhu, J. Mol. Liq., 2024, 409, 125499.
24. J. Li, C. Ren, K. S. Moon and C. P. Wong, 2019 IEEE 69th Electronic Components and Technology Conference (ECTC), 2019, pp. 2296–2301.
25. D. Xie, Q. Pang, Y. Zhao, Y. Li, F. Li and H. He, J. Appl. Polym. Sci., 2022, 139, 52940.
26. D. Guzmán, X. Ramis, X. Fernández-Francos and A. Serra, Polymers, 2015, 7, 680–694.
27. T. Okabe, Y. Oya, K. Tanabe, G. Kikugawa and K. Yoshioka, Eur. Polym. J., 2016, 80, 78–88.
28. Y. Oya, M. Nakazawa, K. Shirasu, Y. Hino, K. Inuyama, G. Kikugawa, J. Li, R. Kuwahara, N. Kishimoto, H. Waizumi, M. Nishikawa, A. Waas, N. Odagiri, A. Koyanagi, M. Salviato and T. Okabe, Chem. Phys. Lett., 2021, 762, 138104.
29. S. Maeda and K. Ohno, J. Phys. Chem. A, 2005, 109, 5742–5753.
30. S. Maeda, K. Ohno and K. Morokuma, Phys. Chem. Chem. Phys., 2013, 15, 3683–3701.
31. K. Ohno and S. Maeda, Chem. Phys. Lett., 2004, 384, 277–282.
32. K. Ohno and S. Maeda, Phys. Scr., 2008, 78, 058122.
33. T. Okabe, T. Takehara, K. Inose, N. Hirano, M. Nishikawa and T. Uehara, Polymer, 2013, 54, 4660–4668.
34. Y. Kawagoe and T. Okabe, Comput. Mater. Sci., 2023, 228, 112333.
35. T. Watanabe, Y. Kawagoe, K. Shirasu and T. Okabe, Int. J. Solids Struct., 2023, 283, 112489.
36. N. Odagiri, K. Shirasu, Y. Kawagoe, G. Kikugawa, Y. Oya, N. Kishimoto, F. S. Ohuchi and T. Okabe, J. Appl. Polym. Sci., 2021, 138, 50542.
37. Y. Kinugawa, Y. Kawagoe, K. Oine, K. Ryuzono, Y. Hoshikawa and T. Okabe, Adv. Composite Mater., 2024, 33, 1420–1432.
38. Y. Kinugawa, Y. Kawagoe, K. Shirasu, K. Ryuzono and T. Okabe, Int. J. Mech. Sci., 2025, 306, 110779.
39. T. Watanabe, Y. Kawagoe, Y. Hoshikawa, Y. Nakai, K. Ryuzono and T. Okabe, Int. J. Solids Struct., 2025, 308, 113130.
40. Y. Kawagoe, K. Kawai, Y. Kumagai, K. Shirasu, G. Kikugawa and T. Okabe, Mech. Mater., 2022, 170, 104332.
41. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
42. S. L. Mayo and B. D. Olafson, J. Phys. Chem., 1990, 94, 8897.
43. A. K. Rappe and W. A. Goddard III, J. Phys. Chem., 1991, 95, 3358–3363.
44. L. Zhang, S. V. Zybin, A. C. V. Duin, S. Dasgupta, W. A. Goddard and E. M. Kober, J. Phys. Chem. A, 2009, 113, 10619–10640.
45. B. Demir and T. R. Walsh, Soft Matter, 2016, 12, 2453–2464.
46. R. W. Hockney and J. W. Eastwood, Computer simulation using particles, CRC Press, 2021.
47. M. Tuckerman, B. J. Berne and G. J. Martyna, J. Chem. Phys., 1992, 97, 1990–2001.
48. W. Shinoda, M. Shiga and M. Mikami, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 134103.
49. M. S. Gordon, J. M. Mullin, S. R. Pruitt, L. B. Roskop, L. V. Slipchenko and J. A. Boatz, J. Phys. Chem. B, 2009, 113, 9646–9663.
50. M. A. Downey and L. T. Drzal, Polymer, 2014, 55, 6658–6663.
51. K.-F. Lin, W.-Y. Shu and T.-L. Wey, Polymer, 1993, 34, 277–288.
52. T. Vialon, H. Sun, G. J. Formon, P. Galanopoulo, C. Guibert, F. Averseng, M. N. Rager, A. Percot, Y. Guillaneuf, N. J. V. Zee and R. NicolaÃ, J. Am. Chem. Soc., 2024, 146, 2673–2684.
53. H. Su, S. A. Jansen, T. Schnitzer, E. Weyandt, A. T. Rösch, J. Liu, G. Vantomme and E. W. Meijer, J. Am. Chem. Soc., 2021, 143, 17128–17135.
54. F. Ugozzoli and C. Massera, CrystEngComm, 2005, 7, 121–128.

---