Heat transfer and SET voltage in filamentary RRAM devices

TL;DR

This paper investigates heat transfer mechanisms in filamentary RRAM devices using COMSOL simulations and analytical models, revealing improved approximations and proposing a universal expression for SET voltage related to material properties.

Contribution

The study introduces refined analytical models for heat transfer in filamentary RRAM, compares them with detailed simulations, and proposes a universal SET voltage expression based on Lorenz number considerations.

Findings

Parallel heat transfer mechanisms improve temperature estimation.

Electrode proximity significantly affects heat transport modeling.

Low Lorenz number materials may enable lower SET voltages.

Abstract

We study the heat transport in filamentary RRAM nano-sized devices by comparing the accurate results of COMSOL modeling with simplified analytical models for two complementary mechanisms: one neglecting the radial heat transfer from the filament to the insulating host, while the other describing the radial transport through the dielectric in the absence of the filament heat transfer. For the former, we find that the earlier assumed simplification of the electrodes being ideal heat conductors is insufficient; a more adequate approximation is derived where the heat transport is determined by the adjacent proximities of the filament tips in the electrodes.We find that both complementary mechanisms overestimate the maximum temperature yet offering acceptable results. However, the two in parallel provide a better analytical approximation. In addition, we show that the Wiedemann-Franz-Lorenz law helps the analysis when the Lorenz parameter is chosen from the actual data. We present an approximate expression for the SET voltage possessing a high degree of universality and predicting that filament materials with low Lorenz numbers can be good candidates for the future low set voltage devices.

Figures

Tables

Table 1. TABLE I: Material parameters used in numerical modeling [ 1 , 2 , 3 , 4 , 5 , 6 ]

Material	TiN	HfO2	HfO2-x	Hf
Thermal Conductivity[W/K.m]	11.9	0.5	20a	23
Electrical Conductivity[S/m]	106	10-2	105	3$\cdot {10}^6$
Specific Heatc [J/kg.K]	545.33	120	130a	144
Relative Permittivity	-$\infty$b	25	-$\infty$a ,b	-$\infty$b
Density[kg/m3]	5220	9680	12000a	13310

Table 2. TABLE II: Device dimensions used in numerical modeling

Device/Dimension (nm)	$H$	$h$	$l$	$d$
Device I	30	10	100	6
Device II	100	50	100	20

Equations

$$-\nabla (\kappa (T)\nabla T)=\mathbf{J}\cdot \mathbf{E}.$$

$${\nabla }^{2}T=-\frac{\sigma }{\kappa }{E}^2\equiv -\frac{E^2}{L}\frac{1}{T}.$$

$$T(z)=T_j+\frac{E^2{h}^2}{8L{T}_j}\left[1-{\left(\frac{2z}{h}\right)}^2\right]$$

$$\delta T=\frac{E^2{h}^2}{8L{T}_j}$$

$$T_r=(T_j-T_0)\frac{d}{2r}+T_0$$

$$T_j=T_0+\frac{1}{2}\sqrt{T_0^2+\frac{{\kappa }_f}{{\kappa }_m}\frac{E^{2}hd}{4L}}$$

$$T_{\max }=T_j+\delta T$$

$$T(z)=T_{\max }\exp \left[-{\left({\mathrm{erf}}^{-1}\left(\frac{\mp z}{T_{\max }}\sqrt{\frac{2E^2}{\pi L}}\right)\right)}^2\right]$$

$$\frac{T_{\max }}{T_j}\mathrm{erf}\left(\sqrt{\ln \left(\frac{T_{\max }}{T_j}\right)}\right)=\mp \frac{h}{2T_j}\sqrt{\frac{2E^2}{\pi L}}.$$

$$\frac{1}{\rho }\frac{d}{d\rho }\rho \frac{d\overline{T}}{d\rho }+{\beta }^2(T_0-\overline{T})=0\mathrm{with}\beta =\sqrt{2\frac{{\kappa }_m}{{\kappa }_{i}Hh}}$$

$$\overline{T}(\rho )=\frac{V^{2}d}{4h{\kappa }_i\beta (h/{\sigma }_f+2H/{\sigma }_m)K_1(\beta d/2)}{K}_0(\beta \rho )+T_0$$

$$R_{th1D}$$

$$R_{\mathrm{thrad}}$$

$$V_{\mathrm{SET}}=h\sqrt{\frac{\delta \mu }{{\sigma }_f{\tau }_T}}.$$

$$V_{\mathrm{SET}}\approx \sqrt{\frac{\kappa \delta \mu }{\sigma C}}.$$

$$V_{\mathrm{SET}}\approx \sqrt{\frac{L}{3}}T.$$

$${\int }_0^{(T^{\prime }(z){)}^2}d{\left(\frac{dT}{dz}\right)}^2=-\frac{2E^2}{L}{\int }_{T_{\max }}^T\frac{dT}{T}.$$

$$\frac{dT}{dz}=\pm \sqrt{\frac{2E^2}{L}}\sqrt{\ln \left(\frac{T_{\max }}{T}\right)}$$

$${\int }_1^{T/T_{\max }}\frac{dt}{\sqrt{\ln (1/t)}}=\pm \frac{1}{T_{\max }}.\sqrt{\frac{2E^2}{L}}{\int }_0^{z}dz.$$

$$-2{\int }_0^{\sqrt{\ln (T_{\max }/T)}}\exp (-y^2)dy=\pm \frac{1}{T_{\max }}.\sqrt{\frac{2E^2}{L}}z.$$

$$T(z)=\frac{T_{\max }}{\exp \left({\left({\mathrm{erf}}^{-1}\left(\frac{\mp z}{T_{\max }}\sqrt{\frac{2E^2}{\pi L}}\right)\right)}^2\right)}.$$

$$\frac{T_{\max }}{T_j}\mathrm{erf}\left(\sqrt{\ln \left(\frac{T_{\max }}{T_j}\right)}\right)=\mp \frac{h}{2T_j}\sqrt{\frac{2E^2}{\pi L}}.$$

$$-{\kappa }_i\left[\frac{1}{\rho }\frac{\partial }{\partial \rho }\rho \frac{\partial }{\partial \rho }T(\rho ,z)+\frac{{\partial }^2}{\partial z^2}T(\rho ,z)\right]=0.$$

$$\frac{h}{\rho }\frac{d}{d\rho }\rho \frac{d\overline{T(\rho )}}{d\rho }+\left[\frac{\partial T(\rho ,z)}{\partial z}{|}_{z=h/2}-\frac{\partial T(\rho ,z)}{\partial z}{|}_{z=-h/2}\right]=0.$$

$$-{\kappa }_i\frac{\partial T(\rho ,z)}{\partial z}=-{\kappa }_m\frac{\partial T(\rho ,z)}{\partial z}.$$

$$\frac{h}{\rho }\frac{d}{d\rho }\rho \frac{d\overline{T}(\rho )}{d\rho }+2\frac{{\kappa }_m}{{\kappa }_i}\left(\frac{T_0-\overline{T}(\rho )}{H}\right)=0.$$

$${\rho }^2\frac{d^2\overline{T(\rho )}}{d{\rho }^2}+\rho \frac{d\overline{T}(\rho )}{d\rho }-{\beta }^2{\rho }^2(\overline{T}(\rho )-T_0)=0.$$

$$\overline{T}(\rho )=C_1{I}_0(\beta \rho )+C_2{K}_0[\beta \rho ]+T_0.$$

$$\overline{T(\rho )}=C_2{K}_0(\beta \rho )+T_0.$$

$$-{\kappa }_f\frac{\partial T(\rho ,z)}{\partial z}=-{\kappa }_i\frac{\partial T(\rho ,z)}{\partial z}.$$

$$C_2=\frac{-1}{\beta K_1(\beta \rho )}\frac{d\overline{T}(\rho )}{d\rho }\mathrm{and},$$

Full text

I Introduction

Modern resistive random access memory (RRAM) devices are about 100-300 nm thick including two metal electrodes and the the insulator layer of about 10-70nm thickness. For the the insulator thickness of $\sim 10$ nm and and bias voltage of $\sim 1$ V, the electric field strength in the insulator $E\sim 10^{8}$ V/m. In response to that field, a conductive path through the insulator can be created in the form of a conductive filament. The regime of continuous conductive filament corresponds to the device ON state created at the SET transition and discontinued at the RESET transition. When the current flows through the filament, a significant Joule heat is liberated creating a temperature distribution that affects the device operations. Here, we present an analysis of that temperature distribution with the emphasis on its underlying physics facilitating the device understanding and design. This is achieved by (1) numerically generating the temperature distribution by means of COMSOL multiphysics package and (2) comparing the latter with two analytically solved limiting cases: (a) heat transport to the electrodes dominated by the metal filament without lateral spreading , vs. (b) that dominated by the lateral spreading of thermal energy preceding its dissipation at the electrodes. Based on our conclusions, we present a related expression for the filamentary RRAM SET voltage that possesses a high degree of universality.

II Numerical Modeling

We use an axillary symmetric model of a RRAM device, the cross-section of which is shown in Fig.1; a description of its COMSOL modeling routine is presented in Appendix A. The dielectric layer corresponds to the metal oxide HfO2. (As widely accepted, the oxygen sweeps creating oxygen vacancies that form a conductive filament of a non-stoichiometric composition HfO2-x.) Following a number of implementations [3, 7], we choose TiN for the metal electrode materials. The material properties are listed in Table I.

Two device dimensions chosen for modeling are described in Table II. The top electrode is connected to a source voltage of 0.5 V while the bottom one is grounded. An ambient temperature of 300K is maintained at the top and bottom. The side walls of the device are taken to be electrically and thermally insulated, which reflects an array geometry with zero inter-device currents. With those boundary conditions, COMSOL solves the coupled heat and electrostatic equations producing the temperature distributions presented in Fig. 2(d).

As shown in Figs. 2(a), 2(b), 2(c), and 2(d), the temperature distributions are a maximum of about 610 K and 576 K at the center of the filament for the Devices I and II respectively. The temperature decay scales in the axial and lateral directions are comparable, although the functional forms are different. The similarity of the two decay scales is due to the mutually balancing competing factors: a better thermal conductivity in a metal vs. the greater area facing the dielectric. In conceivable cases of different material parameters or geometrical dimensions, one of those factors can dominate.

III Limiting cases

For comparison, we present below the two limiting cases corresponding to (a) the 1D heat transfer along the filament and (b) the dominating radial transport of heat from the filament to the insulator. At the end, we combine the two limiting cases as the thermal resistance in parallel to obtain a single heat transport model.

III-A 1D Thermal Transport

Assuming the thermal conductivity of the insulating layer much less than that of the filament material, the heat will flow only through the filament corresponding to a 1D problem. Consider a filament of length $h$ embedded in an insulator sandwiched between two metal electrodes and subjected to an external voltage V as shown in Fig. (1). The ends of the filament are at the filament-junction temperature $T_{j}$ which depends on the ambient temperature $T_{0}$ maintained at the surface of the metal electrodes.

Consider a steady state heat equation with Joule heat as a source term,

[TABLE]

where, $\bf{J}$ is the current density and $\bf{E}$ is the electric field in the filament. Since $\delta\kappa/\kappa\ll\delta T/T$, we treat thermal conductivity to be a constant; a broader analysis is presented in [8]. Using the Ohms law, $\bf{J}=\sigma\bf{E}$, the above equation takes the form,

[TABLE]

Here we have accounted for the Wiedemann-Franz-Lorenz Law[9] ,[10], stating that $\sigma/\kappa=1/LT$ where $L$ is the Lorenz number. Sommerfeld[11], showed that for ‘good’ metals, $L=(\pi k_{B})^{2}/3e^{2}=2.44\cdot 10^{-8}$ W$\Omega$K*-2* where, $k_{B}$ is the Boltzmann’s constant and $e$ is the electron charge.

In reality, the measured $L$ varies between different metals metals and can be temperature dependent, off from the Sommerfeld’s value within an order of magnitude. [12] In particular, for the case under consideration, the values from Table I lead to $L\approx 6.67\cdot 10^{-7}$ W$\Omega$K*-2*, significantly different from the Sommerfeld’s prediction. That may indicate that the filament material is not a homogeneous metal representing perhaps a mixture of metallic and insulating phases [13]. Earlier, the Wiedemann-Franz-Lorenz relation was used in modeling of phase change memory devices [14].

Eq. (2) has a simple solution for the regime of low heat, $\delta T\ll T$, predicting the parabolic temperature distribution,

[TABLE]

with the maximum (center) temperature by,

[TABLE]

above the junction (at $z=\pm h/2$) temperature $T_{j}$.

The latter solution was earlier presented in [8] and [15]. It was assumed [15], that the junction filament temperature is the same as the ambient temperature, $T_{j}=T_{0}$, which may seem intuitively justified because of the high thermal conductivity of the metal electrode . An important correction here is that $T_{j}$ can be significantly higher than $T_{0}$ because the thermal transport in the electrode is dominated by the geometrically small region adjacent to the filament tip that plays the role of the point heat source at the electrode interface. Indeed, the thermal resistance of the series of elemental co-centric semi-spherical layers centered at the tip, is determined by the small distance contributions, $r\sim d$. The significant difference between $T_{j}$ and $T_{0}$ is illustrated in Fig. 3.

Analytically, the temperature distribution at distance $r$ in the electrode from the above point source satisfying the condition that $T(r)=T_{0}$ at $r\rightarrow\infty$, takes the form,

[TABLE]

The corresponding heat flow through a hemispherical surface of radius $d/2$ in the electrode, $2\pi(d/2)^{2}\kappa_{m}\nabla T$ with $T$ from Eq. (5) must be equal to that supplied by the filament, $(2\delta T/h)(\pi d/2^{2})$ as follows from Eq. (3). That equality yields,

[TABLE]

The maximum temperature is now given by,

[TABLE]

with $T_{j}$ from Eq. (6) and $\delta T$ from Eq. (4).

For numerical estimate, we use $L\approx 6.67\cdot 10^{-7}$ W$\Omega$K*-2*, $T_{0}=300$ K, and $V=0.5$V, yielding $\delta T\equiv T_{\textrm{max}}-T_{j}\approx 91$ K and $94$ K respectively for Device I and Device II, while the junction temperatures are $T_{j}\approx 514$ K and $495$ K respectively, which is close to the numerical solution presented in Figs. 2(a) and 2(b).

As shown in Appendix B and the earlier general treatment [8], Eq. (2) can be integrated more accurately to result in

[TABLE]

where erf*-1* is the inverse error function, and the maximum temperature is defined by the transcendental equation

[TABLE]

These results are plotted in Figs. 2(a) and 2(b) predicting $\delta T\equiv T_{\textrm{max}}-T_{j}\approx 85$ K and 80 K, for Device I and Device II respectively, rather close to the above COMSOL result.

III-B Radial temperature distribution

In spite of the relatively small thermal conductivity of an insulator, the heat flowing from the filament to the metal electrodes through the insulator can be significant because of the large areas of the filament/dielectric and dielectric/electrode interfaces. To understand the corresponding radial distribution of temperature, we adopt a simplified model of a very thin insulating layer where the transversal temperature variation can be neglected, thus approximating $T$ with its average value $\overline{T}(\rho)$.

Equating the radial component of the divergence of heat flow $\kappa_{i}h\rho^{-1}(\partial/\partial\rho)\rho(\partial/\partial\rho)\overline{T}$ to the heat $\kappa_{m}(\overline{T}-T_{0})/H$ absorbed by the electrodes leads to the equation,

[TABLE]

where $\beta$ plays the role of the temperature reciprocal decay length.

The solution to Eq. (10) is expressed through the Bessel function with the coefficient determined by the condition that the heat flux $\pi hd(d\overline{T}/d\rho)$ at $\rho=d/2$ equals the filament generated Joule heat. As shown in Appendix (C), this results in

[TABLE]

where $K_{0}$ and $K_{1}$ are the modified Bessel functions of order 1 and 0.

For numerical estimate, we use the above presented device parameters, which yields the temperature decay length $1/\beta=2.51$ nm and 10.24 nm for Device I and Device II respectively. The resulted temperature distribution inside the insulator is plotted in Fig. 2(c) along with the numerically solved (COMSOL) temperature distribution. A somewhat overestimated temperature of the filament is due to the approximation neglecting thermal exchange through the filament bases.

III-C Parallel Transport

A more accurate temperature estimate can be obtained by considering the parallel connection of the thermal resistance from 1D heat transport $R_{th1D}$, and that of radial heat transfer, $R_{thrad}$, which are derived from respectively Eq.(4) and Eq.(11),

[TABLE]

where $R$ is the electrical resistance. The parallel connection results in a smaller total thermal resistance, $1/R_{\textrm{th}}=1/R_{\textrm{th1D}}+1/R_{\textrm{thrad}}$.

The maximum temperature can be calculated using $\delta T=R_{\textrm{th}}\cdot V^{2}/R$. For numerical estimate of Device I and Device II, we use the above presented parameters, which yields $\delta T\equiv T_{\textrm{max}}-T_{j}\approx 82$ K and 79 K, respectively. In Sec. III, we obtained $T_{j}$ to be 514 K and 495 K for Device I and II respectively, using which the maximum temperature, $T_{\textrm{max}}$, of the filament is found to be 596 K and 574 K which is yet closer to the COMSOL numerical solution.

IV SET Voltage

Our thermodynamic analysis of resistive switching relates $V_{\textrm{SET}}$ as in [[16], eq.(8)],

[TABLE]

where $\delta\mu$ is the difference in chemical potentials (per volume) between the insulating phase and the conducting phase of the filament, and $\tau_{T}$ is the thermalization time. In the approximation of 1D heat transfer, one gets $\tau_{T}=h^{2}/\chi$. The thermal diffusivity, $\chi=\kappa/C$ where $C$ is the volumetric heat capacity. Hence $V_{\textrm{SET}}$ can be written as,

[TABLE]

Using the Einstein model, $C=3nk_{B}$ where $n$ is the number density ($N$/$Vol$) of atoms. The chemical potential can be expressed in terms of chemical potential per atom ($\delta\mu_{a}$) as $\delta\mu=n\delta\mu_{a}$. Estimating $\delta\mu_{a}\approx k_{B}T$ where $T$ is the temperature at which the filament phase is formed, and adopting the notations of Widemann-Franz-Lorenz law the SET voltage is then expressed as,

[TABLE]

For numerical estimate, we take the above discussed $T=600$K and $L=6.67\cdot 10^{-7}$ W$\Omega$K*-2* which yields $V_{\textrm{SET}}\sim 0.3$ V which is in the ballpark of measured SET voltages [7, 17]. Perhaps more importantly, Eq. (16) predicts that RRAM devices in which the filament materials have low Lorenz numbers can operate at the correspondingly lower SET voltages.

V Conclusion

We have presented the accurate COMSOL modeling of thermal transport in filamentary devices in comparison with simplified analytical models, one neglecting the radial heat propagation, while the other approximating the temperature distribution as uniform in transversal direction. Our results show that for the typical device parameters, both approximations somewhat overestimate the filament temperature, although the errors ($\lesssim$ 10%) are not very significant. In particular, we have shown that the model of 1D (along the filament) thermal transport is surprisingly accurate when amplified with a realistic heat transport through the electrode suggested here. These simplified models can serve as a convenient express analysis tool. Our predicted values of filament temperature fall in the ballpark of the earlier measured and modeled values [18, 19, 20]; some differences can be attributed to the particular structural and material parameter choices, and the effects of interfacial resistances [14] neglected here.

Also, we have demonstrated that the Wiedemann-Franz-Lorenz law can be used in device analysis when the value of Lorenz parameter is taken to correspond the experimental data on electric and thermal conductivity, which can result in significant deviations from its Sommerfeld’s value.

Finally, we have presented an approximate expression for the SET voltage in filamentary RRAM structures, possessing high degree of universality and correctly predicting the measured values. That expression points at the filament materials wilt low Lorenz numbers as candidates for the future low switching voltage devices.

VI Acknowledgement

We are grateful to I. V. Karpov, R. Kotlyar, and V. I. Kozub for useful discussions.

Appendix A COMSOL Model

Our COMSOL algorithm is as follows.

1. Open the Model Wizard.

2. Choose 2D Axisymmetric as **Space Dimension

**3) Choose 2D AC/DC module and add Electric Currents sub module in Physics.

4. Choose Heat Transfer Module and add Heat Transfer in Solids submodule in Physics.

5. Create the Geometry of the MIM structure as in Fig(1).

6. Create Blank Materials in the Materials node and add material parameters given in Table I to create the required materials.

7. Assign the materials to the corresponding domains.

8. The Heat Transfer in Solids submodule will have four different necessary default subnodes - add Temperature boundary condition and select the top boundary of the top electrode and bottom boundary of the bottom electrode and choose 300K in the user defined temperature section.

9. The Electric Currents submodule also has four different necessary default subnodes - add Electric Potential boundary condition, select the top boundary of the top electrode and set 0.5 V in the Electric Potential box, add Ground boundary condition, select the bottom boundary of the bottom electrode.

10. In Multiphysics node select all the domain in Electromagnetic Heat Source sub-node to couple the Electric Currents submodule and Heat Transfer in Solids submodule,

11. Create Mesh.

12. Select Study.

13. Obtain results in desired form from the Results node.

Appendix B 1D Heat Transfer

In this appendix, we present the solution to the 1D heat transfer problem discussed in Section III-A. Multiplying Eq.(2) with $dT/dz$ and integrating yields,

[TABLE]

or,

[TABLE]

where the origin is at the center of filament. Define $t\equiv T/T_{\textrm{max}}$. Integrating Eq. 18 from $t(0)=1$ to $t(z)=T/T_{\textrm{max}}$ gives the temperature distribution in the filament.

[TABLE]

With a substitution of $y^{2}=\ln(1/t)$, the integration changes to,

[TABLE]

The LHS of the Eq. (20) is an error function. Hence we get the following temperature distribution after the integration.

[TABLE]

The maximum temperature can be calculated by integrating Eq. (18) from $t(0)=1$ to $t(h/2)=T_{j}/T_{\textrm{max}}$. The integration results in a transcendental equation of the form,

[TABLE]

which can be solved numerically.

Appendix C Radial Heat Transfer

In this appendix, we present more in detail derivation of Eq. (10) for the radial heat transfer and its solution in Section III-B. Writing the heat equation in cylindrical coordinate,

[TABLE]

Since, $h/H\ll 1$ the ratio $\delta T/T\ll 1$ along $z$-direction, hence replacing $T(\rho,z)$ by average temperature $\overline{T(\rho)}$ and integrating along z-axis from $z=-h/2$ to $z=h/2$, Eq.(23)

[TABLE]

The second term of the LHS of the Eq.(24) can be reduced to averages after applying the Neumann boundary condition at the metal-insulator interface. The condition states that,

[TABLE]

The temperature gradient in the electrode is replaced by the slope $\Delta T/\Delta z$ and also by using the average temperature $\overline{T}(\rho)$ as the boundary temperature of the metal-insulator junction, Eq.(24) becomes,

[TABLE]

Eq.(26) can be rewritten as,

[TABLE]

The solution to Eq. (27) is the linear combination of zeroth order modified bessel functions $I_{0}$ and $K_{0}$,

[TABLE]

where, the constants $C_{1}$ and $C_{2}$ is determined by the Neumann boundary condition. We omit $I_{0}$($C_{1}=0$) from the solution as it increases with $\rho$. Hence the radial temperature distribution is given by,

[TABLE]

The Neumann boundary condition becomes,

[TABLE]

The filament receives power $P=V^{2}/R$ from the external bias source, where the electrical resistance $R$ is given by sum of electrodes and filament resistances $R=h/\sigma_{f}\pi r_{f}^{2}+2H/\sigma_{m}\pi r_{f}^{2}$, assuming the current takes the shortest path through electrode to filament having the same cross section area as the filament. The power generated by the filament is given by the product of heat flux ($-\kappa_{f}dT/d\rho$) and the surface area of the cylinder at $r=d/2$. Using the boundary condition and equating the power received and generated we obtain,

[TABLE]

Substituting the latter equations, the radial temperature distribution in the insulator becomes,

[TABLE]