Effects of quantum recoil forces in resistive switching in memristors

1. Introduction. Memristors—electric devices whose resistance “remembers” the history of the previously passed current Chua_Memristor —have attracted a lot of theoretical and experimental attention over the past decade, notably since the first working memristive device was manufactured in 2008 Memristor2008 . Memristors are particularly appealing in the context of their possible use for in-memory computing and hardware neural networks, as opposed to the von Neumann architecture Ielmini_inmemory . Current-controlled resistive switching in such devices typically occurs due to ion migration in a dielectric layer between the anode and the cathode resulting in the growth or rupture of conductive filaments Wang_memristors_review . Interestingly, despite considerable progress made on the technological side, there are still gaps in the understanding of the physical processes underlying resistive switching and retention of the resistive states in memristors. The reason is the stochasticity of the filament dynamics and possible interplay of various scales and effects, such as Joule heating, surface diffusion, strong electric fields, etc. Lee_review ; Ielmini_NatComm2019 An improved theory or simulation frameworks are particularly desirable when one descends to the nanoscale: the constrictions, i.e., the thinnest segments of the conductive filaments in the corresponding devices, allow for ballistic (quantum) charge transport and exhibit quantized conductance at room temperature, with the quantization being an attractive property to encode discrete states of the memristor Xue_review2019 ; Gao_Gquant ; Minnekhanov_OrgElec2019 ; Shvetsov_GquantPPX ; Kharlanov2022 . At the same time, interaction of quantum currents with the classical filament and its molecular environment constitutes a challenging task to describe, limiting rational design of such memristors.

In this letter, by no means aiming at a complete description of the dynamics of quantized-conductance filaments, we analyze a possible role of electron-induced, “recoil” forces acting on the ions Dundas_NNano2009 ; Todorov_PRB2010 in these dynamics, which is also potentially related to the extended lifetimes of filaments with well-quantized conductance values. The point is that conductance quantization in quasi-1D systems in the units of $G_{0}=e_{0}^{2}/\pi\hbar\approx 77.5~{}\mu\text{S}$ ($e_{0}$ is the elementary charge, $\hbar$ is the Planck’s constant) is not quite fundamental, being typical of “smooth”, adiabatic filaments Landauer ; Buttiker ; in principle, it does not have to hold for irregular filaments forming in memristors. However, the back reaction of quantum electrons passing through the filament could push it toward smoother shapes consistent with minimum recoil, integer electron transmission probability (transmittance), and thus quantized conductance. Indeed, in our recent study Kharlanov2022 , we identified that recoil forces should be large enough at voltages $V\sim 1~{}\text{V}$ and conductances near $G_{0}$ to compete with the interatomic forces, potentially affecting the evolution of the filament. Together with that, our experiments showed that the thinnest filaments with the conductance $G=G_{0}$ and even $G_{0}/2$, where the quantum recoil forces should be most pronounced, demonstrated the longest renention, i.e., stability with time Kharlanov2022 . However, Ref. Kharlanov2022 described conduction electrons using a semi-quantitative continuous free-electron model and lacked a dynamical simulation of the filament stochastic evolution, aiming at an order-of-magnitude estimation of the effects. Moreover, interatomic interactions in the filament were simply reduced to the latter having a certain amount of surface energy, like a liquid drop. In fact, since ion migration in memristors, e.g., surface diffusion, proceeds as a sequence of hops Ielmini_NatComm2019 , the lack of the dynamics strongly limits the model used in Ref. Kharlanov2022 , especially when it comes to simulate an actual resistive switching event. In contrast, in the present work we develop an atomistic quantum theory of electron recoil forces and incorporate it into a molecular-dynamics (MD) framework with realistic interatomic potentials to observe a quantum analogue of electromigration Electromigration triggering a resistive switching. The theory and simulations, described in the following sections, demonstrate that if the conditions for the quantized conductance are not satisfied, resistive switching can be facilitated by the growing recoil forces. This new resistive switching mechanism, which is important for memristors with atomic-scale filaments and high current densities, represents the main novelty of the present work, providing an alternative to the conventional mechanism due to ion migration in the electrostatic field between the cathode and the anode Lee_review ; Ielmini_NatComm2019 ; Xue_review2019 . Our results thus indicate that quantum recoil should be taken into account when describing memristive nanodevices and could potentially be used to achieve better retention characteristics of their quantized-conductance states.

2. The model. Similarly to earlier studies Dundas_NNano2009 ; Todorov_PRB2010 of current-induced forces, we use a semiclassical, Born–Oppenheimer (adiabatic) description of the nanoconstriction supporting a quantum current: the positions of the ions are treated classically, while the conduction electrons are described within the framework of nonequilibrium Green’s functions (NEGFs) applied to the tight-binding model Do_NEGFreview . Namely, the Hamiltonian of the system reads:

where $\hbar$ and $m$ are the Planck constant and the electron mass, $\mathbf{X}=(\mathbf{X}_{1},\ldots,\mathbf{X}_{N})$ and $\mathbf{P}_{i}=M_{i}\dot{\mathbf{X}}_{i}$ are the coordinates and the momenta of the ions, respectively, $U(\mathbf{X})$ is the “bare” ion-ion potential, and $u(\mathbf{x};\mathbf{X})$ is the effective (Kohn–Sham) potential. For simplicity, we assume one orbital per atom, so that the electron field operator with spin projection $\sigma=\pm 1/2$ is expanded as $\hat{\psi}_{\sigma}(\mathbf{x})\approx\sum_{i}\hat{c}_{i\sigma}e_{i}(\mathbf{x}-\mathbf{X}_{i})$ over an atomic-orbital basis set $\{e_{i}(\mathbf{x})\}$, with the coefficients $\hat{c}_{i\sigma}$ being the electron annihilation operators. The transfer integrals and on-site energies for these orbitals form the $t_{ij}$ matrix. The Ehrenfest equations for the atomic coordinates give $\dot{\mathbf{P}}_{i}\approx-\mathbf{\nabla}_{i}U(\mathbf{X})+\mathbf{F}^{(\text{el})}_{i}(\mathbf{X})$, where the force acting on atom $i$ due to electrons equals:

with the coefficients $\mathbf{f}^{i}_{jk}=-\int e_{j}^{\ast}(\mathbf{x})\mathbf{\nabla}_{i}u(\mathbf{x};\mathbf{X})e_{k}(\mathbf{x})\mathrm{d}^{3}x$. Note that the total forces acting between the atoms contain the term coming from conduction electrons both in the presence and in the absence of the current through the system. The coefficients $\mathbf{f}^{i}_{jk}(\mathbf{X})$, albeit not expressible in terms of $t_{ij}(\mathbf{X})$, can be estimated under reasonable assumptions. Namely, let us assume that (i) the Kohn–Sham potential can be represented as a sum of even contributions $u_{i}(\mathbf{x}-\mathbf{X}_{i})$ from individual atoms, (ii) the basis orbitals, in turn, possess a definite parity, and (iii) $u_{i}$ decays with a characteristic spatial scale $1/\zeta$ between atom $i$ and its neighbors. Then $\nabla_{i}u(\mathbf{x};\mathbf{X})\sim\zeta u_{i}(\mathbf{x}-\mathbf{X}_{i})\cdot\frac{\mathbf{x}-\mathbf{X}_{i}}{|\mathbf{x}-\mathbf{X}_{i}|}$, and $\mathbf{f}^{i}_{ii}$ vanishes identically due to parity. The nonnegligible coefficients that remain are $\mathbf{f}^{i}_{jj}=-\nabla_{i}t_{jj}(\mathbf{X})\sim\mathbf{f}^{i}_{ij}=\mathbf{f}^{i\ast}_{ji}\sim-\zeta\mathbf{n}_{ij}t_{ij}$, where $\mathbf{n}_{ij}\equiv\frac{\mathbf{X}_{j}-\mathbf{X}_{i}}{|\mathbf{X}_{j}-\mathbf{X}_{i}|}$, for the nearest neighbors $j$ of site $i$ (denoted below as $j\in\text{NN}(i)$). The $\zeta$ parameter is naturally expected to be around the inverse atomic radius. The resulting rough estimation of the force thus reads:

The expectation values in Eqs. (3), (4) are taken over a nonequilibrium many-body state describing a current of conduction electrons. In order to describe it using NEGFs, we introduce a customary partition of the atomic sites into two “leads” $\text{L},\,\text{R}$ and the filament F connecting them:

so that $\hat{c}_{\text{L},\text{R},\text{F}}$ denote columns with all annihilation operators in the leads or the filament, $H_{\text{L},\text{R},\text{F}}$ are the Hamiltonian matrices of these three subsystems, and $V_{\text{L},\text{R}}$ couple the filament to the two leads. The retarded NEGF of the filament is introduced in a standard way,

where the self-energies $\Sigma_{\text{L},\text{R}}=V_{\text{L},\text{R}}^{\dagger}g_{\text{L},\text{R}}V_{\text{L},\text{R}}$ stem from isolated-lead Green’s functions $g_{\text{L},\text{R}}(E)=(E-H_{\text{L},\text{R}}+\mathrm{i}\eta)^{-1}$, and $\eta\to+0$. The transmission probability in the form of Caroli is Caroli :

with the broadening operators $\Gamma_{\text{L},\text{R}}=\mathrm{i}(\Sigma_{\text{L},\text{R}}-\Sigma_{\text{L},\text{R}}^{\dagger})=2\pi V_{\text{L},\text{R}}^{\dagger}\rho_{\text{L},\text{R}}V_{\text{L},\text{R}}$, where $\rho_{\text{L},\text{R}}(E)$ are the spectral functions of the isolated leads.

Now, to find the expectation value in Eq. (3), we introduce “incident” wave functions $\varphi^{(n)}_{\text{L},\text{R}}$ in the isolated left or right leads, with energies $E_{n}$, so that the scattering state $\psi^{(n)}$ of the whole system is given by NEGFs: e.g., for the waves coming from the left lead, $\psi^{(n)}\equiv\bigl{(}\psi_{\text{F}}^{(n)},\;\psi_{\text{L}}^{(n)},\;\psi_{\text{R}}^{(n)}\bigr{)}=\bigl{(}G_{\text{F}}V_{\text{L}}^{\dagger}\varphi^{(n)}_{\text{L}},\;(\mathds{1}+g_{\text{L}}V_{\text{L}}G_{\text{F}}V_{\text{L}}^{\dagger})\varphi^{(n)}_{\text{L}},\;g_{\text{R}}V_{\text{R}}G_{\text{F}}V_{\text{L}}^{\dagger}\varphi^{(n)}_{\text{L}}\bigr{)}$. The rest of the calculation is straightforward and takes its simplest form if both atoms belong to the filament, $j,k\in\text{F}$:

where $f(E-\mu)$ is the Fermi–Dirac distribution in the lead with chemical potential $\mu$. Observing that the term in braces gives nothing but the spectral function $\rho_{S}(E)$, we readily yield:

Used together with expression (3), this yields the electron-induced, “recoil” force:

where $\mu_{\text{L}}=E_{\text{F}}-(1-\beta)e_{0}V$ and $\mu_{\text{R}}=E_{\text{F}}+\beta e_{0}V$ are the chemical potentials of the leads, $V$ is the applied bias, $E_{0}$ and $E_{\text{F}}$ are the conduction band bottom and the Fermi level, respectively. The $\beta$ parameter controls a possible asymmetric voltage drop in the leads attached to the constriction Xue_review2019 ; in what follows, we will assume the symmetric case $\beta=1/2$ unless otherwise specified. Note that the energy integral in Eq. (10) extends to the whole conduction band, unlike the case of the Landauer–Büttiker formula for the conductance Landauer ; Buttiker :

However, as we noted above, the zero-current interatomic interaction (e.g., the force field used in MD simulations) should include both the bare term $-\mathbf{\nabla}_{i}U$ and the one given by Eq. (10) integrated over $[E_{0},\;E_{\text{F}}]$. Therefore, the extra, “renormalized” force $\mathbf{F}^{(\text{el})}_{i}(V)-\mathbf{F}^{(\text{el})}_{i}(0)$ acting on the atoms at a nonzero bias originates from the electronic states with $E\in[\mu_{\text{L}},\;\mu_{\text{R}}]$, similarly to the Landauer formula.

It should also be noted that the mere momentum conservation law (which can be derived from the field-theoretic expression (3)) gives the total force acting on the constriction aligned with a certain axis $z$:

with the sum running through open channels in the leads with directions $\hat{z}_{\text{L},\text{R}}=\mp 1$, wave numbers $k_{n}^{\text{L},\text{R}}(E)$, and transmittances $\mathcal{T}_{n}(E)$. This force vanishes for a perfectly ballistic conductance ($\mathcal{T}_{n}=1$); moreover, for identical leads, the integral cancels out except the $[\mu_{\text{L}},\;\mu_{\text{R}}]$ segment, analogously to the Landauer formula:

For a single open channel at the Fermi surface and a small bias, $F^{(\text{el})}_{z}\approx\frac{2}{\pi}e_{0}Vk_{\text{F}}(1-G/G_{0})$, where $k_{\text{F}}$ is the Fermi momentum, which suggests that forces measured in tenths of $\text{eV}/\text{\AA}$ can act on the thinnest filaments with their transmittances not close to unity.

Intuitively, it appears that similarly to the total force (13), the recoil forces (10) acting on individual atoms should also be suppressed in the ballistic regime with the transmittance $\mathcal{T}(E)$ close to an integer, which, according to Eq. (11), is characterized with a quantized conductance. To test this physical intuition, below we consider an analytically solvable case of a 1D chain of atoms, followed by an MD simulation of a nanofilament with realistic interatomic interactions and the quantum recoil on top of them. While the former system is meant to serve as an idealized “toy model”, letting us estimate the magnitude of the forces, their dependence on the bias, and correlation with the transmittance, the latter simulation will give a definitive answer on how the recoil forces can affect the evolution of a conducting filament and drive a resistive switching in a memristor.

3. A quantum point contact in a 1D atomic chain. Let us number the atoms of such a chain, schematically shown in Fig. 1(b), with $j=\pm 1,\pm 2,\ldots$ (right/left leads, respectively) and $j=0$ (the single central atom comprising the “filament”), so that at zero bias, the equilibrium positions of the atoms are proportional to the lattice parameter $a$, $\mathbf{X}_{j}=ja\mathbf{e}_{z}$, along the $z$ axis. The nearest-neighbor transfer integrals involving the central atom are $t_{0,\pm 1}=v_{\text{R},\text{L}}\in\mathbb{R}$, while all the other $t_{j,j\pm 1}=t<0$; let us also take the on-site energies $t_{jj}=\epsilon\delta_{j,0}$. As mentioned, within such a toy model, we are interested in estimating the electron-induced force acting on the central atom, treating all the other, lead atoms as fixed, and the electron transmittance of the chain. Now, working within the NEGF formalism described above, it is straightforward to find the “incident-wave” states, $\varphi_{\text{L},\text{R};j}^{(k)}=\sqrt{2a/\pi}\sin jka$, $E_{k}=2t\cos ka$, $0<k<\pi/a$, and the self-energies $\Sigma_{\text{L},\text{R}}(E)=v_{\text{L},\text{R}}^{2}\int_{0}^{\pi/a}\mathrm{d}{k}\;\phi_{\text{R};1}^{(k)}\phi_{\text{R};1}^{(k)\ast}/(E-E_{k}+\mathrm{i}\eta)=(v_{\text{L},\text{R}}^{2}/2t^{2})(E-\mathrm{i}\sqrt{4t^{2}-E^{2}})$. After that, the filament Green’s function (a $1\times 1$ matrix in our case) reads $G_{\text{F}}(E)=(E-\epsilon-\Sigma_{\text{L}}-\Sigma_{\text{R}}+\mathrm{i}\eta)^{-1}$, together with Eq. (7) giving a closed-form expression for the transmittance:

where $\chi\equiv(v_{\text{L}}^{2}+v_{\text{R}}^{2})/2v_{\text{L}}v_{\text{R}}$. The conductance at small biases, i.e., near the Fermi surface $E\approx E_{\text{F}}=0$, is $G=G_{0}T(E_{\text{F}})=G_{0}/(\epsilon^{2}t^{2}/4v_{\text{L}}^{2}v_{\text{R}}^{2}+\chi^{2})$. Careful evaluation of the total force acting on the central atom ($j=0$) gives an expression:

where $g_{\text{R}}^{1,1}(E)=(E-\mathrm{i}\sqrt{4t^{2}-E^{2}})/2t^{2}$ and the approximation (4) is used. One immediately observes that at zero bias, the force contains a prefactor of $v_{\text{R}}-v_{\text{L}}$ coming from $\Phi^{(p)}$ and thus vanishes for a symmetric chain, $v_{\text{L},\text{R}}=t$. Below we quote the leading order of the zero-bias force in $v_{\text{R}}-v_{\text{L}}\to 0$ (the on-site energy $\epsilon$ does not contribute within the linear order):

where we have made use of the fact that for our chain at half-filling, $E_{0}=2t<0$ and $E_{\text{F}}=0$. At a finite but small bias, assuming that $v_{\text{L},\text{R}}=(1\pm\gamma)t$, $\gamma\to 0$ (which could be a result, e.g., of a deviation of the central atom from the origin along the $z$ axis), one can find the “renormalized” force:

In fact, it is the reflection symmetry mapping the single-atom “filament” onto itself and resulting in $F^{\text{(el)}}_{0,z}\to-F^{\text{(el)}}_{0,z}$, $V\to-V$, $\beta\to 1-\beta$, $v_{\text{L},\text{R}}\to v_{\text{R},\text{L}}$, which forbids the contribution proportional to $V(v_{\text{L}}-v_{\text{R}})$ for a symmetric voltage drop ($\beta=1/2$). In the next section, we will see that for larger filaments with their different ends coupled to different leads, this cancelation does not occur, leading to pronounced finite-bias recoil effects even for symmetric geometries. An estimation with the interatomic distance $a\sim 2.5~{}\text{\AA}$, $\zeta\sim 1/a$, $V\sim 1~{}\text{V}$, $t\sim-2~{}\text{eV}$, and $v_{\text{L}}-v_{\text{R}},\;\epsilon\sim 0.5~{}\text{eV}$ gives the zero-bias force of about $0.7~{}\text{eV}/\text{\AA}$ and the finite-bias contributions around $0.02\text{--}0.03~{}\text{eV}/\text{\AA}$. While the latter figure is quite small, the former one is considerably large at the scale of interatomic interactions and agrees with our above expectations of the total recoil force (13).

It is also interesting to note the effect of the applied bias on the stability of the position $Z_{0}$ of the central atom. Indeed, assuming that the transfer integrals involving this atom are $v_{\text{L},\text{R}}\sim te^{\mp\zeta Z_{0}}$ and $\epsilon=0$, one has the equilibrium position $Z_{0}=0$ with both the electronic force $F^{\text{(el)}}_{0,z}(V)$ and the ionic contribution $-\partial U/\partial{Z_{0}}$ vanishing (due to $v_{\text{L}}=v_{\text{R}}$ and the symmetry, respectively). If the central atom infinitesimally deviates from $Z_{0}=0$, the total force acting on it becomes $F_{0,z}(V;Z_{0})\approx-2\zeta^{2}tZ_{0}\bigl{[}1+8/\pi+(\beta-1/2)e_{0}V/\pi t\bigr{]}-(\partial^{2}U/\partial Z_{0}^{2})Z_{0}$. Now, taking the ion-ion term in the form of a screened Coulomb interaction with the nearest neighbors, $U=q_{\text{eff}}^{2}\bigl{(}|Z_{0}+a|^{-1}+|Z_{0}-a|^{-1}\bigr{)}+(\text{terms independent of }Z_{0})$, where $q_{\text{eff}}$ is the screened ion charge, we arrive at the following hessian:

While we assume that $\partial F_{0,z}/\partial Z_{0}<0$ at zero bias to ensure the stability of the chain before applying the voltage, the latter can make it unstable if $\operatorname{sgn}{V}=\operatorname{sgn}(1/2-\beta)$, above the critical value:

In particular, estimations show that the destabilizing zero-bias term in brackets $1+8/\pi\approx 3.5$ can well compete with the stabilizing Coulomb one $2q_{\text{eff}}^{2}/\zeta^{2}a^{3}|t|\approx 5.7(q_{\text{eff}}/e_{0})^{2}$ for partially screened charges. For example, for $q_{\text{eff}}=0.8e_{0}$ and $\beta=0$, application of a voltage above $V_{\text{crit}}\approx 1.7~{}\text{V}$ results in a localized counterpart of the Peierls instability affecting the central atom’s position Peierls . In fact, this voltage-induced instability could become even more pronounced for more complex filament geometries, in particular, not exhibiting a reflection symmetry.

To complement the discussed analytical estimations, Fig. 1(a) demonstrates the total recoil force in a chain with $\epsilon=0$ versus the conductance (i.e., transmittance), while Fig. 1(b) resolves the force and the charges at different atoms around the quantum point contact. Note that according to Eq. (3), the recoil forces arise from interactions of the ions with Friedel oscillations of the electron density Friedel , thus, these forces should decay into the bulk together with the density oscillations. The amplitude of the forces, as we observe in Fig. 1(b), is large enough to perturb the geometry of the contact and either push it toward a configuration with a transmittance close to unity (when the recoil is negligible), or to facilitate ion migration. In this regard, it is worth noting that even for a constriction $\sim 10~{}\text{\AA}$ (four atoms) long, the conventional, electrostatic force driving ion migration in memrsitors should roughly be $\sim 0.1~{}\text{eV}/\text{\AA}$ at the voltage of $1~{}\text{V}$—compared with $\sim 1~{}\text{eV}/\text{\AA}$ in Fig. 1(b)— which highlights the potentially important role of the recoil forces in the evolution of thinnest conductive filaments in memristive devices. Note in this connection that, although we have mentioned above that the recoil force should be renormalized in some way to avoid double counting upon its addition to the “bare”, classical ion-ion interaction $U(\mathbf{X})$, the oscillations in Fig. 1(b) reflect the quantum nature of the force and cannot be renormalized out after subtraction of the interference-free, classical term. Let us now switch to a more realistic numerical simulation of the effect of the recoil forces, to complement the outcomes of the present paragraph (in particular, to observe the renormalized recoil forces exceeding $0.1~{}\text{eV}/\text{\AA}$).

4. Simulation of a Cu filament with realistic interatomic forces. Of course, the filaments actually forming in memristors are neither one atom wide, as in Fig. 1(a), nor centrosymmetric: they result from stochastic ion migration, e.g., during the so-called electroforming of the memristor, which precedes its normal operation Lee_review . Therefore, the effect of the recoil forces on such, general filaments at room temperature may be beyond the conclusions drawn from the above 1D toy model at $T=0$. To simulate the dynamics and study the recoil-driven ion migration along such filaments, we developed a custom MD code implementing the embedded-atom model (EAM) of the interatomic Cu-Cu forces Sheng_EAM and used a pre-generated typical initial geometry of a conductive filament between two ideal fcc electrodes [Fig. 2(a)], further warmed up to room temperature to obtain a more representative configuration. Following the adiabatic approximation, the recoil forces were calculated for instantaneous ion positions using Eq. (10) and then renormalized by subtracting their zero-bias values. Approximation (4) was used for the force coefficients $\mathbf{f}^{i}_{jk}$ with $\zeta=1/a$, where $a=2.54~{}\text{\AA}$ is the nearest-neighbor distance in the bulk fcc copper. When parameterizing the tight-binding model (2), the differences between the on-site energies of the atoms were neglected, while the transfer integrals were evaluated as $t_{ij}=-1.8~{}\text{eV}\times e^{-\zeta|\mathbf{X}_{i}-\mathbf{X}_{j}|}$, approximating the corresponding DFT results for $4s$ states in Cu dimers. As for the Green’s functions $g_{\text{L},\text{R}}$ of the semi-infinite leads, explicit Fourier-series expressions in the nearest-neighbor approximation were used for them, allowing one to avoid known poor conditioning issues PoorConditioning . Following a typical MD simulation methodology, we first minimized the energy of the system at absolute zero, then warmed it up to $T=300~{}\text{K}$, and finally ran a $150~{}\text{ps}$ production MD, using the Langevin thermostat. Being quite time-consuming, evaluation of the recoil forces was done each $0.1~{}\text{ps}$ of the simulation, while the Langevin equation was integrated with a $1~{}\text{fs}$ timestep.

The recoil forces calculated for the initial geometry [Fig. 2(a)] exhibit a strong asymmetry, despite the atomic positions for this configuration were intentionally chosen at the sites of an ideal fcc crystal. The asymmetry is caused by the bias $V=+1~{}\text{V}$ applied to the anode (the right electrode in the figure): in agreement with Eq. (3), the forces point toward the cathode, which has a higher electron density due to the excess electron waves arriving there from the lead. The forces are measured in several tenths of $\text{eV}/\text{\AA}$, as in the case of the 1D atomic chain. As finite temperature is switched on, the conductance, initially being very close to $3G_{0}$, develops fluctuations around its quasisteady state and further features a sharp resistive switching around $t=110~{}\text{ps}$ [Fig. 2(b)]. During all this time, the total recoil force remains directed toward the cathode, facilitating ion migration in this direction; in agreement with our above analysis, the force drops down when new conduction channels open. The result of the ion migration is seen in the initial and the final states of the filament presented in Fig. 2(b): a whole atomic layer near the anode gets rarified, with the atoms adding to the central part of the filament and thus enhancing the overall conductance. Note that, though the migration has the same direction as that of $\text{Cu}^{+}$ ions driven by the electric field, the latter is absent in our simulations, and the ion transport occurs solely due to momentum transfer from conduction electrons we (perhaps somewhat misleadingly) refer to as the “recoil” here. Thus, we conclude that a specific type of ion electromigration occurring due to a nonequilibrium quantum current of charge carriers can drive resistive switching in the atomic-scale-filaments regime.

5. Discussion. Finally, a couple of words should be said on the applicability of the approximations used above. First, the tight-binding model is a useful, easily parameterizable framework, but is more suitable for describing semiconductors, rather than metals. At the same time, in atomic-thickness filaments, the screening effect may be weaker than that in the bulk metal, with the real physics lying between the tight-binding model and the free-electron one considered recently in Ref. Kharlanov2022 . Second, in the context of these two models, a nontrivial thing to be noted regards the contributions of individual atoms to the total Kohn–Sham potential $u(\mathbf{x};\mathbf{X})$: the additive representation of the latter used in the present work clearly deserves further study and generalization. Third, in the above, we treated electron-current-induced forces in the mean-field fashion: this approximation rests on the fact that at the voltage of $1~{}\text{V}$, around $5\times 10^{14}$ electrons per second pass through a filament with conductance $G_{0}$, so that the time it takes a copper atom to shift by $0.1~{}\text{\AA}$ is enough for around 15, i.e., quite a lot of electrons to pass through the filament. Finally, for transition metals, including copper, localized $d$ orbitals are important, so that to achieve a better numerical accuracy, a multiscale ab initio+MD methodology appears more appropriate. For these further steps, however, the calculations within our simplified model may serve as a proof-of-principle study.

Last but not least, we would like to outline the novelty of the presented study. The main result is the quantum recoil-driven resistive switching mechanism demonstrated above, which works independently of the conventional ion migration in the electrostatic field. Namely, the forces acting on the ions and mediated by electrons get imbalanced due to the difference in the chemical potentials of the two leads, providing an additional driving force for ion migration. Note that unlike the continuous-medium approach applied in Ref. Kharlanov2022 to atomic-thickness filaments, the use of the more realistic atomistic potentials does not lead to the transversal capillary-bridge instability of the filament due to surface tension forces CapillaryInstabilities ; instead, the ion migration proceeds along the filament. Clearly, further simulations with larger and/or asymmetric filaments could reveal more interesting results, however, we leave them beyond the present work because of their higher computational demand. Also, though the very tools we have used for describing the quantum recoil forces—the NEGF formalism and the tight-binding model—are quite conventional and well-established Dundas_NNano2009 ; Todorov_PRB2010 ; Do_NEGFreview ; Caroli , the forces are, to the best of our knowledge, included in a dynamical description of the conducting filament for the first time.

6. Conclusion. To summarize, we have studied the effect of the “recoil” forces between the conduction electrons and the ions comprising a conductive filament on the geometry of the latter. As both analytical estimations and the MD simulation have revealed, the effect of the recoil should be quite pronounced for conductive filaments several atoms wide, like those formed in memristors exhibiting the conductance quantization effect. Unlike the classical volume and surface forces acting on the filament, the quantum recoil is intrinsically related to the transmittance and additionally destabilizes the filaments with the non-integer conductance. Our study thus provides a framework and highlights the importance of taking into account the charge-carrier recoil in the rational design of current and next generations of memristive nanodevices.

O.K. thanks Anton Minnekhanov for fruitful discussions. The work was supported by the Russian Foundation for Basic Research (Project No. 20-07-00696). This work has been carried out using computing resources of the federal collective usage center Complex for Simulation and Data Processing for Megascience Facilities at NRC “Kurchatov Institute”.