## 2 Model and Materials  
### Model  
The geometry of the simulated mushroom PCM cell is shown in Figure 1. The mushroom cells are more relevant for the PCM industry because of the higher density storage and greater susceptibility to thermoelectric heating due to asymmetric geometries. The cell consists of a W (tungsten) top electrode contact, a TiN (titanium nitride) shared top electrode, a GST phase change layer, a TiN narrow bottom electrode (BE), and a W plug surrounded by the thermal and electrical insulation (SiO\({}_{2}\)) with dimensions provided in the Figure 1. The TiN electrode, which has a lower thermal conductivity, acts to confine the energy in the GST layer. An ideal square-wave electrical pulse is biased on the W top electrode contact while the W bottom electrode contact is grounded, which is referred to as the positive bias pulse. A small pulse is applied to calculate the electrical resistance of the cell.  
In this study, we adopt the finite difference method, in which a model space is discretized into three-dimensional Cartesian grids according to the complex geometry of PCM cells with arbitrary cell configuration. The model includes four individual sub-models: electrical, thermal, phase change, and percolation, for a complete description of switching operations in the PCM cells.  
The spatial electrical potential distribution in the cells are modeled by iteratively solving the current continuity equation for each mesh element, and considering particularly the additional current generated due to the Seebeck effect. The calculated electrical current in each mesh element induces local heating as a heat source known as Joule heat. Subsequent heat diffusion changes the temperature distribution, which is calculated by solving the heat conduction equation. Because of the significant contribution compared with Joule heating, thermoelectric heating, consisting of Thomson heat within the GST layer and Peltier heat at the GST-BE interface, should also be considered.  
The current continuity (1) and heat conduction (2) equations are modified to include the thermoelectric effects and are given by  
\[\
abla\cdot J=-\
abla\cdot(\sigma(\
abla V+S\
abla T))=0, \tag{1}\]  
\[dC_{p}\frac{dT}{dt}-\
abla\cdot(\kappa\
abla T)=\frac{J\cdot J}{\sigma}-TJ \cdot\Delta S-TJ\cdot\frac{dS}{dT}\
abla T, \tag{2}\]  
where \(J\) is the current density, \(\sigma\) is the electrical conductivity, \(V\) is the electric potential, \(S\) is Seebeck coefficient, \(T\) is the temperature, \(d\) is the mass density, and \(C_{p}\) is the heat capacity. The thermoelectric contributions are included as thermally driven diffusion current (\(-\sigma\cdot S\
abla T\)) in Equation (1) and thermoelectric heat transport (\(-TJ\cdot\
abla S\), Peltier heating and \(-TJ\cdot(dS/dT)\
abla T\), Thomson heating) in Equation (2) [10].  
The classical theory of nucleation is applied in the phase change model to study switching characteristics according to the calculated temperature distribution [11]. Here, we consider homogenous and heterogeneous nucleation in the bulk material and at interfaces to simulate the crystallization process. The temperature dependent nucleation rate, the growth velocity of the nuclei and the probabilities of crystallization via nucleation or growth are also included. The overall physical property changes of the phase change material, referred to as the percolation effect, are calculated by Bruggeman effective medium approximation [13].