Spatio‐Temporal Correlations in Memristive Crossbar Arrays due to Thermal Effects

Memristive valence change memory (VCM) cells show a strong non‐linearity in the switching kinetics which is induced by a temperature increase. In this respect, thermal crosstalk can be observed in highly integrated crossbar arrays which may impact the resistance state of adjacent devices. Additionally, due to the thermal capacitance, a VCM cell can remain thermally active after a pulse and thus influence the temperature conditions for a possible subsequent pulse. By using a finite element model of a crossbar array, it is shown that spatio‐temporal thermal correlations can occur and are capable of affecting the resistive state of adjacent cells. This new functional behavior can potentially be used for future neuromorphic computing applications.

The ORCID identification number(s) for the author(s) of this article can be found under https://doi.org/10.1002/adfm.202213943.
should not change its resistive state on a scale of at up to 10 years. [4][5][6] Upon the application of appropriate voltage stimuli, at least two resistance states can be established and distinguished, namely the low resistive state (LRS) and the high resistive state (HRS). [4,5] The change of the resistance state is exploited for non-volatile memory applications, for emulating synaptic weights in neuromorphic engineering, as well as for computation-in-memory (CIM). [1,7,8] The latter aspect represents a particular advantage, since such non-Von Neumann technologies not only leads to increased performance, but also means major energy savings. [9,10] In order to achieve very high integration densities, crossbar arrays are used in which the memristive device is sandwiched between two electrodes at their intersection. [10][11][12][13][14] However, this highly scalable architecture has to deal with the problem of sneak-path currents, for which various solutions have already been proposed. [15][16][17] Other parasitic effects that can occur in crossbar arrays are related to thermal phenomena due to Joule heating and become relevant when the switching mechanism of the device is strongly temperature dependent, [4][5][6] which is the case for VCM devices but also for other types of resistive memories. [5,[18][19][20] One spatial effect is thermal crosstalk which arises when the heat of an active device heats up adjacent devices and can thus affect their state. [21][22][23][24][25][26][27][28] Several publications have investigated thermal crosstalk with respect to integration in 2D or 3D arrays. [22,23,26,27] Especially when comparing heat dissipation in different dimensioned arrays, it is found that thermal crosstalk is more challenging in compact 3D arrays as heat accumulates and stays longer in the memristive array due to higher heat capacities. [22,23] In parameter studies, the impact of different sized and shaped filaments, feature sizes, and distances between the memristive cells were analyzed. [21,22,24,25,27,28] Furthermore, the thermal conductivities of the processed materials also influences the heat dissipation in the array. [21,24] The relevance of thermal crosstalk depends on the particular programming case. For example, the resistance state of a half-selected cell in the HRS is much more affected if the more surrounding cells in the LRS are active and thus contribute to higher temperatures. [23] To address the problem of thermal management in memristive arrays, different architectures are investigated and proposed to suppress thermal crosstalk. [23,26,27] In general, one should pay attention on thermal crosstalk when considering high-density memristive arrays. [21][22][23][24]26,27] It should be mentioned that it is quite complicated to approach this thermal

Introduction
The requirements for future electronic devices for information processing are becoming increasingly challenging, as more and more data are to be processed even faster in even smaller devices while consuming as little energy as possible. The major obstacle is that the advancement of conventional complementary metal-oxide-semiconductor (CMOS) technology as the workhorse of the digital age is slowly reaching its physical limits, which results in the end of Moore's law. [1] To overcome this, novel architectures based on memristive devices are considered as promising candidates for future computing solutions and have therefore attracted great attention. [2,3] According to the underlying physical mechanism, one type of these devices is the very well-studied valence change memories (VCMs), which belong to the class of the redox-based resistive random access memory (ReRAM). [4,5] Here, resistive switching occurs due to the migration of oxygen ions resulting in highly non-linear switching kinetics. This is required to solve the voltage-time dilemma where, on the one hand, the memristive device should switch within nanoseconds at a few volts, and on the other hand, a read voltage of a few hundred mV issue experimentally, as it is difficult or even impossible to measure the temperatures and heat fluxes directly. Therefore, a theoretical view, a simulation of finite element, or compact models is required.
In addition, one can observe a temporal effect which occurs due to the transient thermal response of the heated device if one considers time scales similar to the characteristic heating times. [22,23,25,27,29] Thus, with the application of pulse stimuli, the temperature evolution in an active device depends on the previous pulses. [25] Hence, besides the conductance of a memristive cell, the temperature can be considered as a second state variable. As a result, this opens up completely new possibilities for neuromorphic computing. Such short-term temperature dynamics can be used to modulate the resistive state in dependence of the activity. [29] These spatial and temporal effects always occur together and may be of interest for certain applications under defined spatial and temporal conditions. On the one hand, the intention is to suppress these parasitic effects, for example, in terms of security aspects when using a memristive array as classical memory, where unwanted bit flips may occur. [30,31] On the other hand, these effects can be exploited for neuromorphic computing applications. [29] To our knowledge, there is no work yet that deals with spatio-temporal thermal correlations in memristive crossbar arrays in a systematic analysis.
In this work, a finite element model of a crossbar array is presented to investigate time-dependent thermal effects in and between memristive devices. For this purpose, the switching dynamics of a verified compact model of a VCM cell are coupled to the model. With electrothermal simulations, the origin and dynamics of thermal crosstalk and of the thermal accumulation effect are revealed. Furthermore, by using pulse trains, the impact of thermal crosstalk and the thermal accumulation effect with regard to the resistance change is studied. These two thermal effects which occur together under particular conditions then result in the so called spatio-temporal thermal correlations.S1

Methodology
We used the simulation model originally introduced by Staudigl et al. [30] to investigate thermal crosstalk in memristive crossbar arrays. Its limitations are the static approach, thus neglecting the temperature dynamics in transient simulations. In our simulations, the crossbar model is also considered time-dependent, and a dynamic switching model is introduced. The modeling and simulations were performed by the finite element software package COMSOL Multiphysics. [32] Furthermore, we briefly introduce a general method to quantify thermal crosstalk independent of the crossbar structure or the power application.

Modeling of a Memristive Crossbar Array
To analyze thermal effects, a passive crossbar array with one memristive device at each node is chosen since it is probably the simplest structure but also has the best scaling prospects. [12] However, any other memristive device or to address sneak paths even a 1D1R [23,33] or 1S1R [33,34] structure could be used at this point. Figure 1a illustrates the structure of a 5×5 crossbar array laying on a Si/SiO 2 substrate. At each intersection of a bottom electrode (BE) and a top electrode (TE), a memristive device is sandwiched between them, as shown in Figure 1b. As already mentioned, a filamentary VCM device is considered. The filament is formed during an electroforming step [4,35] in a highly insulating metal oxide (MO) layer. According to the temperature-sensitive compact model JART VCM v1b, [36] the memristive device is modeled by a filament consisting of a highly conducting plug region and a disc region serving as the potential barrier where the switching takes place. Due to its relatively high work function, the metal facing the disc region is referred to as the electronically active electrode (AE). The plug region is connected to the ohmic electrode (OE). By applying a negative write voltage, the oxygen vacancies migrate from the OE to the AE. During this SET operation, the device switches from the HRS into the LRS. Due to a reversed voltage polarity, the oxygen vacancies leave the disc region. This leads to the RESET operation and the cell is back again in the HRS. [4,6,36] Thus, the oxygen vacancy concentration is increased in the LRS and decreased in the HRS.
By solving the transient heat transfer equation along with the current continuity equation for the temperature T and the potential φ, respectively, the temperature distribution in the crossbar array and especially in the VCM cell can be determined. Here, k denotes the thermal conductivity, ρ m is the mass density, C p is the heat capacity, σ is the electrical conductivity, E = −∇φ describes the electric field, and j denotes the current density. As thermal boundary condition, the top surface of the crossbar structure is thermally insulating. This assumption was made because simulations with an array encapsulated with, for example, SiO 2 have shown that the differences in temperature are not significant and there is no impact on thermal crosstalk. The assumption is that the bottom side of the substrate, and the vertical boundaries of the model at the end of the electrodes are connected to the ambient temperature T 0 = 293 K representing ideal heat sinks. Thus, the temperatures could even be underestimated, since in a realized chip many layers connected with vias are on top of each other and consequently further away from the heat sinks.
To address the memristive devices, every electrical contact is connected to a certain potential. Upon the application of a V/2 scheme to the crossbar array, a memristive device can be selected. This is referred to as the standard selection pattern. As shown in Figure 1a, only the cell in the center is active by applying V = V app at the third bottom electrode and connecting the third top electrode to ground. It should be noted that only the front faces of the lines are considered as electrical contacts. Besides these electrical contacts, all surfaces are electrically insulating, and for the normal current flow follows j n = 0 A m −2 . Thus, the crossbar structure is electrically not symmetric. The actual crossbar array maintains a certain distance from the model boundary, referred to as padding, to suppress unrealistic effects due to a too close heat sink. Although asymmetry in the work functions and thus in the materials of the electrodes is required for a VCM device, [4] the electrodes are both assumed to consist of platinum to ensure that the thermal crosstalk investigations are not affected by different thermal conductivities of the lines. The geometric and material parameters are shown in Tables 1 and 2, respectively.
To obtain a switching mechanism, the resistive state of the VCM cell in the simulation model is varied by the electrical conductivity of the filament. Following the derivations from Bengel et al., [36] the electrical conductivity in the disc and plug region depends on the oxygen vacancy concentration N disc/plug and can be determined by where e is the elementary charge, VO z is the oxygen vacancy charge number, and µ n is the electron mobility. N disc represents the state variable of the VCM cell whereas N plug is set to a fixed value. Since all other quantities are constant in Equation (3), the oxygen vacancy concentration is a measure for the conductance. The time-dependence of the state variable is modeled as a differential equation and depends on the ionic current I ion : Adv. Funct. Mater. 2023, 33, 2213943  Here, A denotes the area of the filament and l disc the length of the disc region. The further derivation can be found in the publication of Bengel et al. [36] At this point, it is important to note that the ionic current is a function of the vacancy concentration N disc , the voltage across the disc region V disc (for the SET operation), and the maximum temperature of the disc region T disc . The temporal evolution of the ionic current is the origin of the strong non-linearity in the switching kinetics. In contrast to the compact model JART VCM v1b, [36] the temperature can be extracted directly from the finite element model and therefore provides a more accurate description of the reality than the calculation using a thermal resistor. [36] Hereby, the maximum temperature in the disc region is of concern as this is where the switching takes place. Another choice of temperature, that is, the average value, is also conceivable. Note that Equation (4) is solved along with Equations (1) and (2).
In the context of this simulation model, it is referred to a LRS when N disc = 2 × 10 27 m −3 and to HRS when N disc = 1 × 10 24 m −3 . However, during a SET or RESET operation, values in between are also possible. [4] Moreover, by using pulse trains with sub-threshold voltages, a gradual modulation of the resistance is feasible. [37,38] Furthermore, it is also assumed that the RC-times are small enough and thus the charging and discharging of the devices is not considered. This simplification is reasonable, since charging times of down to tens of picoseconds have been observed for such small volumes. [39,40] To focus on the qualitative description of thermal effects, only the SET operation is regarded.

Determination of a Quantity to Describe Thermal Crosstalk
To describe the thermal crosstalk in a memristive crossbar array, time-independent electrothermal simulations of the model are considered in a first step. By using the standard selection pattern, here, the cell in the center is set to the LRS and all others to the HRS. In this manner, the influence of a single cell can be studied. Treating a HRS cell would lead to much less thermal crosstalk because the current density, and thus the Joule heating, is low compared to a LRS cell. As a result, a continuous temperature distribution is achieved, which is shown in Figure 1c. Second, the temperatures in each device are extracted and mapped in a temperature matrix, where T 33 = T sel denotes the temperature of the selected cell. In the third step, a voltage sweep is performed on the simulation model. For each voltage, and thus power dissipation step P i , a corresponding temperature matrix is obtained. It should be noted that only the power dissipation in the disc region of the VCM cell C33 is considered. In the final step, this set of values can be used to determine the thermal resistance R th and heat coupling coefficients α ij . According to Fourier's law, the temperature increase is linearly proportional to the power dissipation. [41] Thus, a linear regression by using and yields to the thermal resistance R th for the selected cell and to the alpha values α ij for all adjacent cells (cf. step four in Figure 1c). Here, P sel denotes the power dissipation in the selected VCM device, and T ij is the temperature in the adjacent cells. From these two equations follows that the heat transfer coefficient is the ratio between the temperature of the adjacent and heated device minus the ambient temperature, respectively: Overall, the alpha values measures the thermal crosstalk independent of the power supply and are characteristic for specified geometries and materials. If no thermal crosstalk can be observed, this results in α ij = 0. The other extreme, the adjacent cell adopts the same temperature as the heated device yields to α ij = 1. Such heat coupling coefficients can be used to investigate thermal crosstalk itself or for the estimation of temperatures in circuit simulations. [30] Likewise, this method is not limited to crossbars and can therefore be applied to other structures as well. For a final evaluation of the influence of thermal crosstalk on the switching behavior, however, one should always also consider the absolute temperatures.

Results and Discussion
In the following, the origin of the thermal crosstalk in crossbar arrays is discussed and a quantitative method to characterize it is presented. Considering transient simulations, the impact of thermally active devices on the resistance state of adjacent cells is illustrated. Furthermore, the thermal accumulation effect observed in dense pulse sequences is analyzed both from a theoretical point of view and in terms of switching dynamics in a simulation model. Finally, spatio-temporal thermal correlations where these two effects occur simultaneously are systematically studied and clarified in an exemplary scenario.

Thermal Crosstalk
For the investigation of the origin of thermal crosstalk, a static electrothermal simulation of a 5×5 crossbar array is considered. According to the standard selection pattern, the LRS cell in the center is selected. Figure 2a   electrodes, with the heat distribution starting from the heated cell C33. Correspondingly, the hot spot in the temperature matrix in Figure 2b is also in the center. The heat generated in the active cell spreads into the attached electrodes and propagates along them. Since the electrodes are only connected to the highly insulating substrate, MO, or to the insulating thermal boundary condition, there is hardly any loss in the heat flux between two devices. At the crossing with another line, the electrodes are only separated by a very thin MO layer containing the filament. For cells C32 and C23 in Figure 2a, heat can flow through the cells into the perpendicular electrodes due to this low thermal insulation. In this way, the heat generated in C33 reaches the diagonal adjacent cell C22. Moreover, to get an idea of the large temperature gradients, Figure 2c illustrates the temperature distribution along the selected bottom electrode. The extremely high temperatures of the hotspot are limited to the small volume of the filament. Figure 2d shows the heat flux along the cross section of the BE. It is clear that the heat flux in the very well conducting electrodes is higher than in lower thermally conducting regions.
To determine the portion of the heat flux that goes into the adjacent devices, the vertical heat flux is monitored along the upper electrode surface. In the figure, it is indicated as a red line. The resulting data in Figure 2e clearly reveal the heat flux into the VCM cell. Reaching the next adjacent cells, the heat flows directly into the MO and thus into the TE. Thereby, the heat flow decreases with the length of the intersection. This will lead to an asymmetric heating in the memristive device (i.e., the right side is hotter than the left side) and therefore could affect the switching. Furthermore, it can be extracted that there must be another source of heat when comparing the left with the right side. Since the current path is only on the left side of the line, a contribution from the electrodes can be identified. Depending on the line resistance, the electrode additionally emits heat due to Joule heating. This also explains the increased device temperatures along the current-carrying part of the selected electrodes (cf. Figure 2b). It should be mentioned that the lower the resistance ratio between the HRS and the LRS is, the larger becomes the effect of self-heating in halfselected cells and thus leads to a third heat contribution. This aspect should be regarded in particular for analog applications.
These findings explain the distribution of the alpha values (cf. Figure 3a) which are derived from Equation (7). In Figure 3b, the alpha values of three different adjacent cells are considered for spacings between 10 and 1000 nm. For all alpha values, an exponential decrease can be observed. A particular characteristic is that the alpha value of cell C32 tends toward a non-zero value for long distances. This is attributed to the Joule heating of the current-carrying part of the selected lines.
In general, thermal crosstalk should be considered for very  small distances, since this effect disappears after a few hundred nanometers. This work will not discuss further parameters, but it is straightforward to imagine that such alpha values can also be derived for the variation of other geometric and material quantities.
In addition, the spacing also has an influence on the thermal resistance of the memristive device and thus on its temperature, which is decisive for the switching kinetics (cf. Figure 3c). This results from the fact that the thermal resistance, as defined here, depends on the thermal properties of its environment. Thus, a small crossbar array structure with densely packaged devices more effectively dissipates the heat to the surrounding VCM cells.
Although steady-state simulations provide important insights into the phenomena, it is also necessary to address the switching dynamics to discover whether the observed temperature increases are only a theoretical influence or whether the resistance state is actually changing. To keep the simulation time on a feasible time scale, transient simulations are performed in a 1 × 3 crossbar array (cf. Figure 3d). Furthermore, this structure is similar to a 1T1R configuration as the memristive devices in a 1T1R configuration are directly connected via exactly one line. The selected cell in the center is in the LRS; the adjacent cells are in the HRS and half-selected. The applied voltage is assumed to be V app = −1.5 V. Increasing or decreasing the voltage would lead to higher or lower temperatures and thus to an acceleration or slowing down of the switching speed. [4][5][6]36] The ten pulses examined have a pulse length of t on = 100 ns and a pause between them of t off = 100 ns (cf. Figure 3e). The latter is long enough to cool down to ambient temperature after each pulse, since the thermal time constant is τ th = 4.5 ns. This will be discussed in detail in the next section. However, this value is extracted directly from the temperature profile of one pulse that reaches the steady state. The rise time respectively fall time is assumed to be t r = t f = 1 ns. The simulation is executed for three different spacings.
In Figure 3e, one can see the transient evolution of the oxygen vacancy concentration and the maximum disc temperature in the two adjacent cells C11 and C13. As expected, the temperatures are higher in the left cell facing the current-carrying line. Moreover, temperatures are decreased for larger spacings between the electrodes. Since the temperature accelerates the VCM switching mechanism, [4][5][6] this leads to the same trends in oxygen vacancy concentration. Thus, the gradual switching process in the adjacent cell depends on thermal crosstalk as well as on the heat contribution of the electrode.
In terms of neuromorphic properties, one can map a longterm plasticity with pulse stimuli on a memristive device. [42] Adv. Funct. Mater. 2023, 33,   Through thermal crosstalk, this capability is additionally transferred to half-selected adjacent cells.

Thermal Accumulation Effect
To better understand the thermal effects over time, it is worth taking a theoretical look at a memristive device subjected to rectangular and unipolar pulse stimulation. During the pulse, the VCM cell heats up while it cools down in the unpowered state. In the most simplified case, the thermal properties of a memristive device could be described through a thermal resistance R th and a thermal capacitance C th . In Figure 4a, the thermal circuit diagram is shown as well as the transient behavior during heating and cooling. In analogy to a RC-circuit, the thermal time constant is defined by where τ th indicates how fast the device heats up or cools down. In case of heating, this is the time after which 63.2 % of the final temperature is reached. [41] According to Newton's law of cooling which can also be derived from Equation (1), the temperature can be expressed as a function of time for heating and cooling Here, T 0 denotes the ambient temperature and T max = T 0 + R th P D the maximum temperature. The latter is determined by R th and the power dissipation P D of the device. Starting from a particular initial temperature T init , the temperature shows an exponential behavior, and therefore, one can distinguish a transient and a steady state. This is often referred to as a transient thermal response. [27] It can be concluded that after 5τ th , the steady state is reached, since the remaining temperature change is less than 1 % of the total change. In Figure 4b, the temperature evolution of a memristive device is illustrated upon the application of two pulses with a pulse width of t on and a pause of t off between them. This temperature profile can be achieved by alternately applying Equations (9) and (10). In the following, ideal voltage stimuli without rise or fall time are assumed and the thermal time constant equals for heating and cooling. It can be recognized that the second pulse is applied before the VCM cell is able to cool down to the ambient temperature. Therefore, the final temperature after the second pulse is higher by ΔT 31 compared to that of the first pulse. When pulses occur in close succession, the temperature in a memristive device can be increased; thus, we refer to this phenomenon as thermal accumulation effect.
Considering an infinite pulse train in Figure 4c with t on = t off , first being in the oscillating transient state, the temperature increases due to the thermal accumulation effect. After the oscillating steady state is reached, the temperature stays between T low and T high and swings around the mean value T mean = (T low + T high )/2. It is obvious that this mean value is also the mean value of T 0 and T max , since pulse width and pause are equal.
By changing the relation between t on and t off in Figure 4d, the oscillating steady state can be shifted to high temperatures when the pulse width is larger than the space. Conversely, it drops for a larger pulse space. In other words, fixing the pulse width will cause the memristive device to get hotter as the frequency is increased. Furthermore, it can be observed that the temperature difference ΔT = T high − T low is larger as there is either more time for heating or cooling.
Since thermal resistance and capacitance are material properties, the thermal time constant can be adjusted, for example, by changing the volume or the switching material. In Figure 4e, the impact of different relations between the thermal time constant and the pulse width (or frequency) are shown. As the thermal time constant becomes large, that is, heating and cooling are slowed down, the oscillating transient state becomes also temporally extended, and in the oscillating steady state, the maxima and minima are closer together. At small τ th , the temperatures oscillate much stronger and less pulses are required to reach the oscillating steady state. To get the same effect with the same device, one has to change the frequency instead of τ th .
As the thermal time constant, in the here studied crossbar array simulated with COMSOL Multiphysics, is under 5 ns, it is expected to see a thermal accumulation effect if the distance between the pulses is lower than 5τ th . Therefore, a single memristive device is investigated with regard to pulse trains with t off ranging from 1 to 20 ns (cf. Figure 5a-c). To measure, how long it takes to reach a certain resistance state regardless of frequency, one can simply count the pulses. Another approach is to define a cumulative pulse time (CPT), which describes the timing in more detail (cf. Figure 5d). According to this, the lengths of the pulses required to switch into a certain resistance state are summed up. If the state is reached within half the pulse length, only half the pulse time contributes to it. Here, the rise and fall time of the pulse is also taken into account.
The transient evolution of the oxygen vacancy concentration and the maximum temperature in the disc region are considered more closely for three different pulse pauses. At first glance, one can see that the temperature peaks start at very low values. Because of the VCM cell being in the HRS at the beginning of the pulse train, in contrast to the theoretical derivation in Figure 4, the power dissipation is much lower and thus leads to smaller possible maximum temperatures. Additionally, there are further deviations in the temperature through the impact of the rise and fall time of the voltage stimuli. These are negligible as long as they are small compared to the pulse length. For t off = 20 ns, one finds that a pulse is hardly influenced by the previous one since t off > 4τ th , and thus, the temperature is almost able to return to the ambient temperature. In this case, it requires 14 pulses to reach N disc = 10 26 m −3 . Due to the thermal accumulation effect, considering the pulse train for t off = 8 ns and t off = 2 ns, only 13 pulses and 11 pulses, respectively, are needed to switch into the same state. Thus, the acceleration in the temperature increase enables faster switching.
The contour plot in Figure 5e summarizes the frequencydependent switching results. The isolines describes which resistance state can be established when applying a certain number of pulses as a function of the pulse spacing. It can be clearly seen that the effect is strongest for small t off . For large pulse pauses, t off > 5τ th , the switching becomes frequency-independent and the thermal accumulation effect can no longer be observed.
Furthermore, it should be emphasized that the thermal time constant of memristive devices depends on many intrinsic and extrinsic factors. It was shown which impact the material properties of the metal oxide have with regard to the thermal time constant. [5] Previous studies also have demonstrated that the thermal time constant varies not only with the array size but also with the type of device that is used at an array node. For example, the thermal time constant of 1D1R device is much higher as that of the standalone 1R cell. [23] With regard to neuromorphic applications, this temporal thermal effect can be used as a second state variable [29] when the distance between two successive pulses is in the order of maximum a few thermal time constants. Thus, the memristive device possibly shows the neuronal ability of a pairedpulse facilitation due to the thermal accumulation effect. Such a short-term plasticity has already been shown for memristive devices in various publications, but has not been attributed to transient thermal behavior. [42][43][44][45] Considering the RESET operation, a reverse voltage is required. Since Joule heating depends only on the magnitude of the potential, paired-pulse depression can also be observed, where switching is accelerated in the other direction.

Spatio-Temporal Thermal Correlations
With the understanding of thermal crosstalk and the thermal accumulation effect, phenomena that occur in crossbar arrays due to different timings can now be addressed. Here, they are referred to as spatio-temporal thermal correlations. An immediate question is how the resistance state of adjacent cells is affected when the pulse spacing is of the order of the thermal time constant. For this purpose, the simulation model of a 1×3 crossbar array from Figure 5d is adapted so that the pulses are modulated by the thermal accumulation effect. Figure 6 shows the chosen parameters and the results when varying the frequency and the electrode spacing. In Figure 6a, the transient behavior of C11 is exemplary given for t off = 1 ns and d = 100 nm. In contrast to the temperature profile of a selected and switching device, the temperature rise is similar to that in the theoretical model (cf. Figure 4). Since this cell is half selected and is initially in the HRS, self-heating is negligible and the heat experienced by this memristive device comes from the selected LRS cell in the center (cf. Figure 6b). Therefore, the switching occurs by the voltage drop of the half-selection, the thermal crosstalk, and is additionally accelerated by the thermal accumulation effect. Figures 6c and 6d indicate the number of pulses required to reach a given resistance state as a function of electrode spacing and pulse pause, respectively, for the adjacent cells. When focusing on a defined pulse pause (one color), the trends are as expected and consistent with the results when looking at thermal crosstalk only (cf. Figure 3e). The difference is that due to the thermal accumulation effect, the adjacent devices change their resistance state faster. It should be noted that the resistance state after the first pulse is independent of the frequency since the thermal accumulation effect only occurs from two pulses on.
To address a more application-related scenario, a 1 × 3 crossbar array is considered next, where individual pulse trains are applied to the top electrodes. Using a ground scheme, only the selected device experiences a voltage drop. Furthermore, the pulse trains are so low frequent that a cell cannot influence itself by the thermal accumulation effect. According to the illustration in Figure 7a, the adjacent cells C11 and C13 are in the LRS and thus are able to generate a large amount of heat due to Joule heating. This affects the temperature in the target cell in the center which is in the HRS. Moreover, the pulses that address the target cell are shifted compared to the adjacent ones. These offsets are chosen to be 2t on , 3t on , and 4t on , which cause the memristive device C12 to be affected by reduced thermal crosstalk. The remaining heat of the adjacent cells enables the target cell to start from a much higher temperature, which is basically nothing more than a temporally increased ambient temperature (cf. Figure 7a-c). Considering the evolution of N disc at a short time offset, the oxygen vacancy concentration increases by three orders of magnitude while it is only about one order of magnitude at a medium time offset and much less than one order of magnitude at a long time offset. This underlines the role of the non-linear switching kinetics in such a memristive device. Due to the thermal crosstalk of the adjacent cells, the temperature is only increased by less than 100 K when the pulse is applied. Here, a minor change of a few Kelvin is sufficient to massively accelerate the switching process.
To conclude, the second state variable of a single resistive device realized by the thermal accumulation effect is, therefore, also accessible in the near environment and not limited to the originating cell. In contrast to the conductance as first state variable, multiple secondary state variables can interact depending on their distance to each other. The combination of thermal crosstalk and the thermal accumulation effect thus leads to spatio-temporal thermal correlations, which can potentially be used for future neuromorphic computing implementations.

Conclusion
In this work, we have presented a finite element model of a memristive crossbar array. In addition, the differential equation of a verified compact model of a VCM cell is coupled with the electrothermal simulation, enabling the simulation of switching dynamics even in large arrays. This simulation model was used to investigate the spatial and temporal thermal Figure 7. a) Simulation setup for a selection scenario in a 1×3 crossbar array, where the pulse train of the target cell C12 is shifted by certain time offset with respect to the adjacent cells. b-d) Transient evolution of the oxygen vacancy concentration N disc and the maximum temperature T disc in the disc region of the target cell C12. The switching is accelerated depending on the given time offset and thus on the remaining heat due to thermal crosstalk. In (c), it is exemplified that the first temperature peak is due to the heat of the adjacent cells and the second peak is caused by the self-heating of the observed device C12. effects that can occur in a highly integrated array. A memristive device is affected by thermal crosstalk across the lines and also by the Joule heating of the electrodes. Furthermore, it has been shown for the first time that spatio-temporal thermal correlations can be observed for device spacings as small as a few hundred nanometers and pulse trains with pauses in the order of the thermal time constant of the memristive device. Based on this effect, novel learning rules can potentially be derived for future neuromorphic computing applications. These findings are likely not limited to crossbar arrays with single VCM devices and can be applied to other temperature-sensitive memristive devices as well, in particular also in 1T1R structures.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.