1. 3D PCM Device Structure and Fabrication

The fabrication process of the vertical 3D PCM device is depicted in Fig. 2. Devices were fabricated on a thermally grown SiO₂/Si substrate. SiO₂ layers were deposited by plasma enhanced chemical vapor deposition (PECVD), while all remaining layers were deposited by physical vapor deposition (PVD) sputtering. A 450 nm base SiO₂ layer was deposited beneath the first bottom electrode (BE) layer to prevent direct contact between the GST/top electrode (TE) TiN and the Si substrate. On top of this base layer, alternating layers of TiN (50 nm) and SiO₂ (100 nm) were deposited and repeated four times to form the stacked BE structure. A Cr hard mask, chosen for its sufficient etch selectivity against both SiO₂ and TiN, was deposited and patterned via photolithography, followed by inductively coupled plasma (ICP) dry etching to define the word line (WL) cut region and electrically isolate adjacent arrays. After Cr removal via ICP dry etching, a 50 nm SiO₂ passivation layer was deposited to electrically isolate the BE sidewalls exposed after the WL cut. This layer prevents unintended electrical contact between the BE layers and the subsequent GST/TiN deposition, ensuring that the current path is confined only through the active hole region. Active hole regions were then defined using the same Cr hard mask approach, followed by ICP dry etching through the full multi-layer stack. After completing the hole etching and Cr hard mask removal, GST for the PCM layer was deposited along the hole sidewalls, and a TE TiN was subsequently deposited. The TE was patterned via ICP dry etching to complete the vertical TiN/GST/TiN cell structure. Staircase BE contacts were formed through sequential photolithography and etching, enabling independent electrical access to each layer.

Fig. 2. Fabrication process flow of the vertical 3D PCM device. (a)–(j) Schematic cross-sections illustrating the fabrication steps. (k) Sequential BE opening from L4 to L1 for independent electrical access to each layer.

The cross-sectional scanning electron microscope (SEM) and transmission electron microscope (TEM) images in Fig. 3 confirm the structural integrity of the fabricated devices across different scales. Fig. 3(a) displays the 1 $\mu\text{m}$ and 5 $\mu\text{m}$ patterns, respectively, confirming that the TiN/SiO₂ multi-layer stack was uniformly patterned regardless of the pattern size. This scalability indicates a stable and reproducible etching process. Furthermore, as shown in the magnified view of Fig. 3(b), the Cr hard mask remained intact after the entire hole etching process, effectively protecting the underlying layers and maintaining a steep sidewall slope of approximately 81°-82°. This confirms that the Cr mask provides excellent etch selectivity, which is crucial for defining vertical channels in a complex multi-layer architecture. The cross-sectional TEM image in Fig. 3(c) illustrates the final vertical architecture of the device. The TEM analysis confirms that the PCM layer (GST) and TE (TiN) layers were conformally deposited along the sloped sidewalls with thicknesses of 25.3 nm and 38.1 nm, respectively. The step coverage of our PVD sputtering system was 40–60% relative to the top surface, which is typical for PVD, posing challenges for precise sidewall thickness control ^[40]. We expect that transitioning to atomic layer deposition (ALD) would address this limitation, providing superior thickness uniformity and more precise control across all vertical layers ^{[30, 41]}.

Fig. 3. Structural characterization of the fabricated vertical 3D PCM device. (a) Cross-sectional SEM images of 1 µm and 5 µm patterns. (b) Magnified SEM image confirming the intact Cr hard mask and well-defined sidewall profile after hole etching. (c) Cross-sectional TEM image confirming conformal deposition of GST and TE TiN.

2. 3D PCM Device Characterization

Electrical measurements of vertical 3D PCM were performed using a LabVIEW program with an 81110A pulse generator, an MSO7104B oscilloscope, a 4155C semiconductor analyzer, and an NI switch system. All electrical measurements were performed on devices with a 5 $\mu\text{m}$ hole diameter. Fig. 4(a) illustrates the I-V characteristics of the vertical 3D PCM for each layer, measured by DC sweep. The devices initially in the high-resistance RESET state (amorphous phase) exhibit an abrupt transition to the low-resistance SET state (crystalline phase) at threshold voltages ($V_{th}$) ranging from 1.5 V to 2.0 V. Even considering the cell-to-cell variation, the devices exhibit stable threshold switching behaviors, which demonstrates the reliable formation of the vertical 3D PCM cells.

Fig. 4. Electrical characteristics of the fabricated vertical 3D PCM device for each layer. (a) I-V characteristics showing threshold switching behavior. (b) R-V characteristics confirming an on/off ratio exceeding two orders of magnitude. (c) SET/RESET switching operation across all four layers under identical pulse conditions (RESET: 5/100/5 ns, SET: 300 ns/5 µs/1.5 µs). Note that SET and RESET states are represented by filled and open symbols, respectively.

As shown in the R-V characteristics of Fig. 4(b), SET and RESET resistances of approximately 8 $\text{k}\Omega$ and 1.3 $\text{M}\Omega$, respectively, were obtained by varying pulse amplitude while keeping pulse timing (rise/width/fall: 5/100/5 ns) constant. This resulted in an on/off ratio exceeding two orders of magnitude, providing a sufficient memory window for multilevel analog weight storage. Furthermore, the switching properties presented in Fig. 4(c) demonstrate that the SET and RESET states exhibit stable switching behavior across all four layers. Notably, the resistance levels extracted from these operations are consistent with the R-V characteristics shown in Fig. 4(b). The consistent results demonstrate that the devices were fabricated with minimal layer-to-layer variation, and that stable device characteristics can be maintained even as the number of layers is increased to achieve higher density. Nevertheless, to further mitigate any remaining cell-to-cell variations, advanced fabrication processes and optimized pulse conditions, such as refined amplitude and width, should be employed for precise device control.

To evaluate the potential for analog synaptic weight implementation, we investigated the multi-state characteristics. As depicted in Fig. 5(a), eight distinct resistance states were successfully achieved and clearly separated. These states were programmed by varying pulse amplitudes, where higher voltage induced greater amorphization of the GST layer, resulting in higher resistance states. These CDF results demonstrate that multi-level resistance states suitable for analog weight representation can be implemented through incremental step pulse programming (ISPP).

Long-term data retention is an important requirement for maintaining trained synaptic weights to ensure reliable weight storage for the inference task. Thus, as shown in Fig. 5(b), we conducted accelerated life testing at elevated temperatures to evaluate the data retention properties. From the resulting Arrhenius plot, all three randomly measured cells demonstrated consistent retention behaviors. By extrapolating this trend, every tested cell can maintain its RESET state for over 10 years at temperatures exceeding 80 °C. This stability ensures that the trained synaptic weights remain reliable over a long time, enabling consistent inference performance in AIMC.

Fig. 5. Multi-level and reliability characteristics of the fabricated vertical 3D PCM device. (a) CDF plots of eight distinct resistance states, with 30 independent measurements per state. (b) Arrhenius plot from accelerated lifetime testing at elevated temperatures (180 °C, 195 °C, and 210 °C). (c) Resistance drift characteristics across four resistance states.

Furthermore, as shown in Fig. 5(c), the resistance drift characteristics of the vertical 3D PCM were investigated. Resistance drift is a well-known phenomenon primarily observed in the amorphous phase, where the resistance continuously increases over time due to structural relaxation ^[21]. The time-dependent resistance shift is typically quantified by the drift coefficient ($\nu$) based on the following power-law relation ^[42]:

where $R_0$ is the resistance measured at time $t_0$, and $\nu$ represents the drift coefficient. As inherent timing delays exist between the programming pulse and the first resistance read-out, $t_0$ was extracted by fitting the resistance–time data to the power-law relation, resulting in $t_0$ values of 0–5 s. In this study, the resistance drift behavior was characterized across four different resistance states. The results show that the drift coefficient for the RESET state is 0.1123, which is consistent with the values reported for conventional GST ^{[17, 21, 22, 43]}. In contrast, the SET state exhibited a significantly lower drift coefficient of 0.0274, demonstrating minimal resistance variation over time. To evaluate the practical impact of these non-idealities, we further investigated how resistance drift affects inference accuracy. The measured drift characteristics were incorporated into a simulation to assess the performance degradation, which will be discussed in the following section.

3. OTS Device Fabrication and Characterization

To verify the characteristics of the OTS device intended for use as a neuron element, a vertical 3D OTS device was fabricated. The same process flow as the vertical 3D PCM device described above was applied, with the only modification being the substitution of GeSe for GST in the deposition step. Electrical characterization was also performed using the same equipment as for the vertical 3D PCM device. The I-V characteristics of the OTS device are shown in Fig. 6(a). The threshold voltage was approximately 3.5 V, measured using a trapezoidal pulse (rise/width/fall: 10/20/10 $\mu\text{s}$). The subthreshold current was measured under DC conditions as shown in Fig. 6(b) (blue curve). The device exhibited a subthreshold current well described by $I = a \cdot \sinh(b \cdot V)$. This compact expression is based on the trap-limited subthreshold conduction model of Ielmini and Zhang ^[44], with $a$ and $b$ treated as effective fitting parameters that encapsulate material- and geometry-dependent factors. A non-linearity coefficient $b = 1.23\text{ V}^{-1}$ was extracted from a global fit to measurements across 8 individual cells, indicating the subthreshold non-linear response required for neuron-element operation. However, this level of non-linearity was found to be insufficient for solving non-linear tasks, as demonstrated in the simulation results presented later, motivating the need to increase $b$. To enhance the non-linearity for neuron applications, an EF process was intentionally applied via a DC forming sweep from 100 fA to 100 $\mu\text{A}$ under a voltage compliance of 1.5 V. This process is known to increase the density of trap states and modify the local amorphous structure of the chalcogenide layer, consistent with the reported EF mechanism ^[45]. Through this process, we obtained a higher non-linearity coefficient of $b = 3.94\text{ V}^{-1}$ from the global fit to post-EF measurements across the same 8 cells, as shown in Fig. 6(b) (red curve). Importantly, threshold switching was confirmed to be preserved after the EF process, as shown in Fig. 6(a), confirming that the EF condition does not induce irreversible changes in the OTS device. These results demonstrate that the 3D OTS device exhibits enhanced subthreshold non-linearity suitable for neuron operation, confirming its potential as a hardware neuron element in the EqProp-based AIMC framework.

Fig. 6. Electrical characterization of the fabricated 3D OTS device. (a) I-V characteristics obtained by trapezoidal pulse (10/20/10 µs), showing the threshold switching before and after EF process. (b) Subthreshold I-V characteristics measured under DC conditions before and after EF process. Scatter plots represent data from 8 individual cells, and solid lines indicate the global fitting curves for after and before EF, respectively.

1. Analog EqProp algorithm

The synaptic weight update in EqProp is determined by the difference in the product of pre- and post-synaptic activities between the free equilibrium (state 0) and the nudged equilibrium (state $\beta$):

where $u_i^0$ and $u_i^\beta$ denote neuron states at the free and nudged equilibria, respectively, $\rho(\cdot)$ is a non-linear activation, and $\beta$ is the nudging strength. This update depends only on local activities, requiring no global gradient signal ^[36].

Kendall et al. ^[37] demonstrated that this mathematical framework can be naturally mapped onto a non-linear resistive network. In this physical implementation, the abstract neuron states ($u$) are represented by node voltages, and the synaptic weights correspond to programmable conductances ($g_{ij}$). Crucially, the "product of activities" in the mathematical rule (Eq. 2) is physically realized through the local power dissipation across the resistors. This allows the weight update to be simplified into a local conductance change driven by the difference in the squared voltage drops between the two operating points (Eq. 3):

where $\Delta V_{ij} = V_i - V_j$ denotes the voltage difference across the synapse.

Building on this AIMC framework, we simulate an analog EqProp framework parameterized based on the electrical characteristics of the fabricated vertical 3D PCM and OTS devices, as described in Section II. PCM synapses are modeled as resistive conductance elements, with the allowable conductance range set by the on-off ratio ($g_{\max}/g_{\min}$) measured from the fabricated device. The fabricated PCM exhibits a conductance window of approximately 0.77 $\mu\text{S}$ to 145 $\mu\text{S}$, and $g_{\max} = 2 \times 10^{-4}\text{ S}$ is adopted as the simulation reference to account for device-to-device variation. For the XOR circuit simulation, the OTS neuron element is modeled using the exponential I-V fit $I = a\sinh(bV)$, where $a$ is the subthreshold amplitude and $b$ is the non-linearity coefficient, extracted from experimental OTS measurements. Although switching-related dynamics such as threshold-voltage variation and temporal effects are not accounted for in this model, the OTS neuron functionality considered in the present EqProp simulation relies on subthreshold non-linearity rather than threshold switching. In the XOR simulation, the OTS bias voltage was selected to position the neuron operation within the subthreshold regime, where the sinh-based compact model reproduces the measured I–V characteristics well. In practical implementations, maintaining a sufficient operating-voltage margin below the threshold-switching voltage would similarly be required to preserve subthreshold device operation.

The analog EqProp framework is validated through two complementary simulation approaches that differ in physical fidelity and computational scale. The small-scale XOR task is solved using Xyce ^[46], an open-source SPICE-compatible circuit simulator, where each neuron node is a circuit node and the free/nudged equilibria are obtained from two separate DC operating-point analyses. For MNIST-scale training, we employ an analytical resistive-network solver that iteratively solves each layer’s equilibrium node voltages via alternating forward–backward sweeps (asynchronous coordinate descent) with fabricated vertical 3D PCM synaptic characteristics ^[47]. To accommodate the non-negative conductance constraint of PCM devices, we adopt the input–output doubling scheme of Kendall et al. ^[37], where each input pixel $p_i$ is expanded to $[p_i, -p_i]$ ($784 \rightarrow 1568$ nodes) and each output class is represented by two nodes, recovering the logit as $\hat{y}_c = z_{2c} - z_{2c+1}$ ($10 \rightarrow 20$ nodes).

2. XOR validation and MNIST Test for Analog EqProp framework

The OTS device characterized in Section II provides the neuron non-linearity required for EqProp-based learning, and the post-EF subthreshold I-V model is used directly in the Xyce circuit simulation of the XOR task. The XOR task requires a non-linear decision boundary and serves as a minimal validation of whether the OTS neuron non-linearity is sufficient to support EqProp-based learning in a two-layer network.

To quantify the role of the non-linearity coefficient $b$ in circuit-level XOR learning, we sweep $b$ from 2 to 5 $\text{V}^{-1}$ and the subthreshold amplitude $a$ under a fixed on-off ratio of 100. For each ($a$, $b$) pair, we evaluate robustness to training hyperparameters by running 27 combinations of learning rate lr, nudging strength $\beta$, and OTS bias voltage $V_{OTS,bias}$ (3 values each), where the XOR success rate is the fraction of these 27 runs that converge to a validation accuracy of 100 % within 30 epochs. The resulting success rates across the ($a$, $b$) parameter space are summarized in Fig. 7(a). The post-EF nonlinearity coefficient ($b = 3.94\text{ V}^{-1}$) of the fabricated device successfully meets the minimum requirement for successful XOR learning under on/off ratio = 100. Fig. 7(b) shows the learning curve for a representative successful run at $b = 4\text{ V}^{-1}$, $a = 4 \times 10^{-6}\text{ A}$, similar to I–V characteristics of the fabricated OTS device ($b = 3.94\text{ V}^{-1}$, $a = 1.66 \times 10^{-6}\text{ A}$). The validation loss, defined as $L = \frac{1}{n} \sum_{i=1}^n (\hat{y}_i - y_i)^2$, decreases monotonically and the accuracy stabilizes at 100 % from epoch 27. These results demonstrate the feasibility of analog EqProp-based XOR learning using the fabricated PCM synaptic and OTS neuron characteristics.

Fig. 7. Effect of OTS non-linearity parameters on XOR learning success at on/off ratio = 100. (a) Success-rate heatmap over the non-linearity coefficient b and subthreshold amplitude a. (27 runs per cell; color encodes the fraction reaching 100 % accuracy). No successes are observed for b ≤ 3 V^-1. (b) Learning curve for fabricated device-similar condition (a = 4×10^-6 A, b = 4 V^-1) at representative hyperparameters (VOTS,bias = 0.3 V, lr= 10^-9 and β = 10^-8), showing validation loss (orange, left axis) and accuracy (blue, right axis), where the shaded region indicates the stable accuracy regime from epoch 27 onward.

While the XOR circuit simulation incorporates both PCM and OTS device characteristics, the MNIST-scale simulation primarily focuses on the PCM synaptic conductance window and resistance-drift characteristics. In this larger-scale analytical simulation, neuron nonlinearity is represented by an ideal diode model. A fully connected MLP with doubled input and output ($784\times 2$ input, 128 hidden neurons, $10\times 2$ output neurons) is trained with the conductance-constrained EqProp model over 10 independent seeds per condition. Fig. 8(a) shows mean test accuracy with standard deviation across on/off ratios from 10 to 1000, with detailed values summarized in Table 1. The dashed line marks the on/off ratio of the fabricated PCM (on/off ratio $\approx 162$). An on/off ratio below 30 causes a substantial accuracy drop of more than 3 percentage points, because an insufficient conductance ratio compresses the effective dynamic range of the synaptic weights. Accuracy saturates above on/off ratio = 100, placing the fabricated device within this regime and confirming that its conductance window is sufficient for analog EqProp training.

3. Inference Accuracy Under PCM Synaptic Non-idealities

To evaluate the influence of programming uncertainty and device-induced conductance variation on inference accuracy, we applied effective multiplicative Gaussian noise to the mapped synaptic weights during inference. This noise represents the combined inference-level effect of programming inaccuracy, device-to-device variation, cycle-to-cycle fluctuation, and conductance variation. For each mapped weight $w$, the noisy weight was defined as $w' = w(1 + \epsilon)$, where $\epsilon$ follows a zero-mean Gaussian distribution with a standard deviation of $\sigma$, so that $\sigma$ corresponds to the relative one-standard-deviation noise amplitude with respect to the target mapped weight. As shown in Fig. 8(b), the mean accuracy drop increases monotonically with the multiplicative noise level, reaching 1.86 percentage points at the largest evaluated noise level. Although the largest noise levels exceed typical relative noise levels reported in PCM-based inference studies ^[24], they confirm the expected trend that larger effective weight/conductance deviations lead to greater accuracy degradation.

Fig. 8. Device-aware simulation of PCM-based analog EqProp on MNIST. (a) Test accuracy versus on/off ratio at fixed $g_{max}$ = 2×10^-4 S. Bars, error bars, and open circles represent the mean, standard deviation, and individual results over 10 independent seeds, respectively. The dashed green line marks the fabricated PCM on/off ratio of approximately 162. (b) Inference accuracy drop under multiplicative Gaussian weight noise during inference, where σ denotes the relative standard deviation of the multiplicative noise term applied to each mapped weight. Error bars indicate ±1 σ across 10 trials. (c) Inference accuracy drop as a function of time under stochastic resistance drift for four variability levels $σ_ν$. Error bars show ±1 σ across 10 trials.

As another important PCM synaptic non-ideality, we further examine post-training inference accuracy degradation under stochastic resistance drift. The state-dependent drift coefficient is assigned to each synaptic conductance based on the nearest measured resistance state in Fig. 5(c), using the experimentally extracted drift coefficients (0.0274-0.1123). Device-to-device drift-coefficient variability is modeled by applying Gaussian dispersion around the assigned state-dependent value, while $\sigma_\nu$ is swept over four levels (0.01, 0.02, 0.05, and 0.10). The simulation is conducted from $t = 1\text{ s}$ to 94,608,000 s (3 years) over 10 trials. Fig. 8(c) shows that the inference accuracy degrades slowly and monotonically in all cases, remaining within 3 percentage points over 3 years across the evaluated drift-variability range. The degradation becomes larger as $\sigma_\nu$ increases, indicating that drift-coefficient variability is an important factor in long-term post-training inference degradation. Within the variability range evaluated in this work, the accuracy degradation remains limited, although stronger drift dispersion beyond the modeled range may lead to additional degradation. The limited degradation observed here may be partly attributed to the energy-minimizing nature of EqProp inference, where inference is formulated as finding the stationary point of an energy function. Under the modeled post-training drift conditions, global conductance shifts can perturb the energy landscape in a partly common-mode manner, allowing the network to converge to a similar equilibrium state. Consequently, the inference process shows limited accuracy degradation under the modeled post-training drift conditions, with a 3-year accuracy drop as low as 0.07 percentage points at $\sigma_\nu = 0.01$.

Notably, these results reveal that accuracy drop is driven primarily by device-to-device variation of the drift coefficient, $\sigma_\nu$, rather than by the drift coefficient itself. When $\sigma_\nu = 0$, the 3-year accuracy drop is negligible regardless of the state-dependent drift magnitude, whereas increasing $\sigma_\nu$ to 0.10 leads to a 2.99 percentage-point drop over the same period.

This identifies reducing $\sigma_\nu$ as the critical mitigation target. Since resistance drift variability tends to scale with the drift coefficient itself, approaches that reduce $\nu$ are expected to also suppress $\sigma_\nu$. Device- and material-level approaches that suppress resistance drift are therefore particularly well-suited to the EqProp framework, which aims to preserve end-to-end analog operation with minimal peripheral overhead. Liner-integrated PCM structures ^{[48- 50]}, doped phase-change materials ^{[51, 52]}, and alternative low-drift phase-change materials ^{[53, 54]} have been shown to reduce the drift coefficient, and represent promising directions for future device optimization. Notably, liner-integrated PCM devices have been reported to reduce not only the drift coefficient but also its device-to-device variability ^{[55, 56]}, directly targeting the $\sigma_\nu$ identified as the primary driver of long-term accuracy degradation. In particular, drift coefficients below 0.02 have been reported ^[22], suggesting that $\sigma_\nu$ could be reduced to levels where long-term accuracy degradation is mitigated.

4. Architectural Comparison with PCM-based AIMC Demonstrations

To clarify the architectural distinction of the proposed approach, Table 2 summarizes representative PCM-based AIMC demonstrations, including recent large-scale PCM computing systems, together with this work. The hardware scale is summarized based on the reported physical PCM device capacity. As shown in Table 2, prior PCM-based AIMC works have demonstrated high task-level accuracy and large physical hardware scales mainly using two-dimensional PCM array architectures combined with mixed-signal or hardware/software execution frameworks that require layer-wise ADC/DAC interfaces. In contrast, this work explores a 3D vertical PCM–OTS framework for fully analog EqProp operation, based on separately fabricated and independently characterized PCM and OTS devices. In this framework, the network state is propagated through analog circuit dynamics without layer-wise ADC/DAC conversion. This distinction highlights the role of the present work as a device-architecture feasibility study toward vertically scalable, high-density, ADC/DAC-free fully analog learning systems, while full PCM–OTS co-integration remains a future validation step.

To quantitatively assess the architectural efficiency of the proposed approach, we estimated both the analog-core energy scale and the intermediate ADC/DAC conversion overhead. Using the mapped PCM conductance scale and characteristic voltage drop ($G \sim 10\ \mu\text{S}$, $\Delta V \sim 1\text{ V}$), the synaptic-array dissipation is estimated to be $\sim 2.0\text{ W}$. Here, $\tau_{\text{inference}}$ denotes the effective integration time for one analog inference, determined by the slowest relevant timescale among network relaxation, OTS response, and output readout. For $\tau_{\text{inference}} = 10\text{--}100\text{ ns}$, this corresponds to an equivalent array-core efficiency on the order of 10 TOPS/W, comparable to the reported analog-tile peak efficiency of 20.0 TOPS/W ^[57]. Further device scaling toward higher-resistance PCM operation, corresponding to lower synaptic conductance, could substantially reduce the synaptic-array power and thereby significantly improve the projected TOPS/W efficiency beyond the present estimate. This estimate is limited to the analog-core level and does not include full chip-level peripheral, routing, or control overhead. Beyond the analog-core estimate, the key architectural advantage of the proposed framework arises from eliminating intermediate ADC/DAC conversions within the internal EqProp relaxation loop. Representative crossbar-based accelerators such as ISAAC ^[58] and PUMA ^[59] identify ADC/DAC interfaces as major peripheral-overhead components, and recent CIM system-level analyses report that ADCs alone can account for approximately 70–80% of total energy and 35–50% of chip area in representative implementations ^{[34, 57, 60]}. These results indicate that signal conversion can be a major bottleneck in conventional mixed-signal AIMC. Therefore, by removing repeated intermediate conversion and signal reloading, the proposed architecture can reduce conversion-related energy, latency, and routing burden as the network depth increases, providing a potential system-level advantage beyond the analog-core efficiency.