Behavioral modeling of integrated phase-change photonic devices for neuromorphic computing applications Open Access

Integrated phase-change photonic devices have recently been shown capable of a quite remarkable range of functionality, including binary and multilevel nonvolatile memory,^1–5 non-von Neumann arithmetic computing,^6,7 and, the focus of this paper, neuromorphic devices and systems.^8–11 Indeed, recently, a small-scale but “complete” all-optical neuromorphic computing system comprising interconnected phase-change photonic neurons and synapses and capable of supervised and unsupervised learning was experimentally demonstrated.⁹ To aid in the future development of larger-scale systems, and so exploit effectively the potential of integrated phase-change photonic devices in the areas of neuromorphic, arithmetic, and in-memory computing, it is important to have access to reliable and accurate behavioral models that can simulate device performance in reasonable time scales. This is because full, physically realistic modeling of integrated phase-change photonic devices and systems is not straightforward, requiring simultaneous handling of interlinked optical, thermal, and phase-switching (amorphization and crystallization) processes, including the effects of any temperature, field, and phase-state induced changes in material/device properties. Usually, such physically realistic approaches to modeling are carried out using finite-element (FE) or finite-difference time-domain (FTDT) methods for the solution of optical and thermal processes, along with nucleation and growth,^12,13 rate-equation,^14,15 or cellular automata approaches¹⁶ for modeling the phase-transformation process itself. However, such physically realistic modeling is computationally resource heavy and the simulation of even the seemingly simplest of processes, e.g., the write/erase process in a binary memory type device, can take many hours using quite powerful computer workstations. Clearly, simulating the performance of even the most basic phase-change photonic computing/processing networks/systems, likely to comprise many hundreds if not thousands of interconnected individual devices, is untenable using such “conventional” methods. A simpler approach is obviously needed, one that encapsulates, in an accurate way, the essential operating characteristics of devices, but one that can be run with modest computing power and produce apparently instantaneous results. It is just such a behavioral model that we outline in this paper.

There are two main types of integrated phase-change photonic device: that based on a straight, rib, or strip type waveguide onto which a phase-change cell is deposited [see Fig. 1(a)] and that based on a microring resonator structure, again onto which a phase-change cell is deposited [see Fig. 1(b)]. The waveguides themselves are often fabricated using silicon nitride, Si₃N₄, on SiO₂ substrates, though the use of the more industry-compatible silicon-on-insulator (SOI) approach (the latter leading to more compact devices and being more inherently CMOS fabrication compatible) is also attractive. Both device types, shown in Fig. 1, can be, and have been, used for binary and multilevel memory,^1–3 and both are suitable for the implementation of photonic neurons^9–11 and synapses.^8–11 A detailed discussion of the operation of these devices for memory and (arithmetic and neuromorphic) processing applications can be found in the literature (e.g., Refs. 1–3, 6–9), and such discussions will therefore not be repeated here. In brief, however, the presence of the phase-change cell on top of the waveguide modifies the optical propagation within the waveguide, and the transmission of light through the waveguide is significantly very different depending on whether the phase-change cell is in its crystalline or amorphous state. For the commonly used phase-change material Ge₂Sb₂Te₅ (or GST for short), for example, the crystalline phase is much more absorbing than the amorphous phase, at the wavelengths (∼1550 nm) typically used in the integrated photonics field. Thus, in the straight waveguide example of Fig. 1(a), the optical transmission through the waveguide is significantly lower when the phase-change cell is crystalline as compared to when it is amorphous, and this forms the basis of the readout process (i.e., a low power optical pulse is sent down the waveguide, and its amplitude on transmission through the waveguide depends on whether the phase-change cell is crystalline or amorphous). To switch the state of the phase-change cell, tailored higher-power (cf. in the readout process) optical pulses are sent down the waveguide to either amorphize (i.e., write to) the cell or to crystallize (i.e., erase) it. The operation of the microring type device shown in Fig. 1(b) is somewhat more complicated than that of the straight waveguide configuration, since changing the state of the phase-change cell on the microring changes the ring’s resonant properties and so changes the amount of light that couples in/out from/to the coupling waveguide.

As pointed out above, any behavioral model for the operation of integrated phase-change photonic devices needs to account for three inter-related processes occurring in the cell, namely, (i) the propagation of light through the waveguide, (ii) the heating of the phase-change cell, and (iii) the crystallization and amorphization (and melting) of the phase-change cell. These three interlinked processes are represented in our behavioral approach via simplified electromagnetic, heat transfer, and phase-change models.

The electromagnetic model uses the effective refractive index of the cell to calculate the power (or amplitude) and phase of the optical propagation mode in the cell, along with the heat generation resulting from optical absorption in the phase-change cell. Since the photonic phase-change devices can be relatively large, with the phase-change cells often several micrometers in length, and since the distribution of crystallized/amorphized regions is not necessarily homogeneous with the cell, we usually split the simulation space into many smaller subcells. The power, ΔP_i, absorbed in the ith phase-change subcell can thus be described as follows:

where P_0i is the power at the leading edge of the subcell, k₀ is the free space wave-vector, L_cell is the length of the subcell, and Im(n_eff)_i represents the imaginary part of the subcell’s effective refractive index. Using this formula, it is possible to calculate the volumetric heat source in each subcell. The optical phase accumulation, Δϕ_i, of the propagating mode on having traversed a subcell is modeled using the following expression:

Using the previous result, it is possible to express the total change in optical phase (ΔΦ_cell) as the difference of the accumulated optical phase for the fully crystalline cell (Φ_cellFC) and the summation of the individual change in optical phase for each subcell (Δϕ_i),

where Re(n_eff)_i is the real part of the effective refractive index of the ith subcell. Reflections at the interfaces between subcells are neglected.

For optical modeling of the microring device of Fig. 1(b), it is also necessary to take account of the transmission and coupling coefficients between the ring and the straight (input-output) waveguide. These are related as follows:¹⁷ 

where b₁ is the mode that passes the ring, b₂ is the mode that enters into the ring, a₁ is the input mode (to the straight coupling waveguide), and a₂ is the mode that couples again to the input/output waveguide after circulating in the ring. The relationship between the mode (b₂) entering into the ring and the mode (a₂) that couples from the ring back into the input/output waveguide is given by

and this enables the derivation of the following expressions for the modes a₂ and b₁ (for a normalized input power of 1 W):

where α and θ are the attenuation and the accumulated optical phase, respectively, after circulating the ring. The coupling and the transmission coefficients, k and t, respectively, are calculated here using approximate analytical expressions.¹⁸ Equations (4)–(7) contain all the necessary information to calculate the amplitude and optical phase of the modes in the ring and input/output coupling waveguide.

The real and imaginary parts of the effective refractive index, n_eff, as used in Eqs. (1) and (2), are accessed by the behavioral model via look-up-tables (LUTs), which, in turn, are obtained by 2D FE eigenmode simulation of the waveguide for various optical powers, temperatures, and phase-state configurations. Such FE simulations essentially provide a library file for a particular device type, whose performance over a vast range of input conditions can subsequently be computed in a fast and efficient manner by the behavioral model. Given the fact that the variation of n_eff is reasonably smooth, we here sampled the temperature dependency between 290 K and 1290 K with a resolution of 100 K, the power dependency between 0 W and 0.02 W with a resolution of 0.002 W, and the crystallinity dependency from 0 to 1 with a resolution of 0.2, and finally, the wavelength dependency was sampled between 1530 and 1560 nm with a resolution of 15 nm. This produced 2178 combinations for the LUT, with intermediate values obtained via linear interpolation (note that the results of the behavioral model did not vary in any significant way if we used higher-order interpolations, inferring that the LUT resolution did not introduce inaccuracies due to under-sampling).

Heat transfer is approached in the behavioral model by using an electrical circuit analogy to solve the heat diffusion (Fourier) equation.¹⁹ Here, the temperature is analogous to voltage, heat to charge, and heat flux to electrical current. We assume that the heat flow is predominantly one-dimensional, in the downward direction perpendicular to the plane of the phase-change cell (i.e., downwards through the ITO, GST, and Si waveguide layers, and into the substrate). This is a reasonable assumption bearing in mind that the ITO and phase-change layers are very thin, compared to their lateral extension, and that their thermal conductivity is much lower than that of the underlying Si waveguide. Some longitudinal heat flow is allowed, however, along the Si waveguide [Fig. 2(a)], reflecting the relatively high thermal conductivity of Si. Thermally, each layer within the device in the vicinity of the phase-change cell is thus modeled as a resistor-capacitor network (RCN), as shown in Fig. 2(b).

The thermal model calculates the temperatures at the top, T_t, bottom, T_b, and middle, T_m, of each layer, with these temperatures being described by the following equations:

where C and R refer to the capacitances shown in Fig. 2(b), with C = ρc_vl and R = l/k, where ρ is the density, c_v is the specific heat, and l is the length (thickness) and k is the thermal conductivity of the relevant material layer. The parameters used in the heat transfer model are provided in Table I.

Turning to the modeling of the phase-change process itself, as pointed out in the Introduction, there are several methods commonly used in the literature that could be thought of as “physically realistic,” in so much as they try to model crystallization (and in some cases, amorphization too) based on physical processes (e.g., related to the competition between volume and surface energies during the formation and growth of a crystal phase). However, these models tend not to be sufficiently efficient computationally for behavioral modeling, in particular, when they have to be interlinked, as in our case, with electromagnetic and thermal simulations. For these reasons, we here use a phenomenological approach based on the well-known Johnson-Mehl-Avrami-Kohnogorov (JMAK) model.^14,20,21 Here, crystalline fraction χ in the cell is given by the following equation:

where K_i is the crystallization rate at the ith time iteration, ν is a frequency factor, E_a is the activation energy for crystallization, k_B is the Boltzmann constant, T_i is the temperature in the ith time iteration, n is the Avrami exponent, and Δt is the time increment between iterations in the model. [Note that although the JMAK approach is used mainly in the literature to model only crystallization, the form of Eq. (11) also allows for the simulation of melting and subsequent amorphization, as explained in Ref. 21.] The crystallization behavior predicted by the JMAK approach is controlled by the parameters E_a [slope in the Arrhenius plot of ln(t_x) vs 1/k_BT, where t_x is the crystallization time] and ν (intersection as T → ∞ in the Arrhenius plot), and so it is obviously important to use appropriate values for such parameters. However, it has been widely observed (see, e.g., Refs. 22 and 23) that a single activation energy cannot adequately describe the crystallization behavior of GST; indeed, in the Arrhenius plot, it is clearly possible to observe two regions that can be fitted using different E_a and ν parameters. For example, Ciocchini et al.²¹ reported E_a values as small as 0.5 eV in the high-temperature region, but much larger values of around 2 eV at low temperatures. Such a change in E_a also affects the position of the intercept of the Arrhenius plot, i.e., the value of the frequency parameter ν. In this work, therefore, we take E_a as being equal to 1.85 eV for T < 400 °C and 0.5 eV above 400 °C. Similarly, ν takes the values 7.35 × 10¹⁹ Hz and 8.1 × 10⁹ Hz below and above 400 °C, respectively. We set the Avrami coefficient, n, equal to 2.6.¹⁴ 

Having described the development of the behavioral model above, we now apply it to simulation of the writing, erasing, and reading of data to/from the phase-change cell. The write process consists of the amorphization of a previously crystalline cell, while the erase process consists of recrystallization of any amorphous region formed during the write process. Readout simply consists of probing, with a suitable low power read pulse, the optical transmission of the device. We begin by considering a straight Si-waveguide device as shown previously in Fig. 1(a) and here with the dimensions as shown in Fig. 3. For simulation purposes, the 4 µm long GST cell is here split into 40 subcells each of 100 nm in length (this resolution was found to be a good compromise between complexity, accuracy, and speed). To write to the GST cell, a relatively high power (here 17 mW) optical pulse of 25 ns duration is sent down the waveguide, as shown in Fig. 4(a) (Multimedia view). The temperature in each of the 40 subcells that make up the GST cell is calculated at each time step (here taken to be 1 ns), along with the resulting fraction of crystallized material in each subcell [see Figs. 4(b) and 4(c) (Multimedia view), respectively]. Finally, the optical transmission, T_r, through the waveguide (equivalent to the readout signal) is calculated, as shown in Fig. 4(d) (Multimedia view). Note that the image shown in Fig. 4 (Multimedia view) is the result of the simulation at t = 45 ns. A link to an animation showing the results generated in the complete simulation can be found in the caption of Fig. 4 (Multimedia view). In each frame of the simulation, we show the temporal evolution of the applied power and the transmittance in Figs. 4(a) and 4(d) (Multimedia view) and a snapshot of the temperature and phase distributions along the cell in Figs. 4(b) and 4(c) (Multimedia view), respectively. From Fig. 4 (Multimedia view), it can be seen that, as expected, the crystalline fraction is reduced on writing (an amorphous mark is written into the cell, but this does not occupy the whole cell volume), and after thermalization, the optical transmission is increased due to the lower absorption (lower k value) of the amorphous phase.

We now turn our attention to the erase process. Here, we use a double-step pulse, i.e., a pulse with an initially high optical power immediately followed by a longer, lower power segment, which has previously been shown to provide good control of the recrystallization process.^3,24 This double-step pulse is applied to a cell that has previously been subject to the write cycle shown in Fig. 4 (Multimedia view). Thus, in Fig. 5 (Multimedia view), we show the result of applying a double-step erase pulse comprising 17 mW for 25 ns, followed immediately by 8.5 mW for 700 ns [Fig. 5(a) (Multimedia view)]. Again, we show the temperature and crystal fraction in each subcell [Figs. 5(b) and 5(c) (Multimedia view)], along with the optical transmission through the waveguide [Fig. 5(d) (Multimedia view)]. By the end of the erase pulse, the cell is, as desired, fully recrystallized (run the associated video file to follow the evolution of crystallized fraction throughout the entire erase process). Note that by varying the amplitude and/or the duration of the second segment of the double-step pulse of the form shown in Fig. 5(a) (Multimedia view), it is possible to set the phase-change cell into various final crystallization states. This gives access to multilevel capabilities needed for neuromorphic applications (e.g., for the setting of synaptic weights), as shown, for example, in Fig. 6 [also note that multilevel states can also be accessed via the degree of amorphization achieved in the write process, controlled by varying the power and duration of the write pulse of the form shown in Fig. 4(a) (Multimedia view)].

As a final example of the capabilities of the behavioral model, we show in Fig. 7 the simulated readout signal for a device having the microring architecture of Fig. 1(b) and with the phase-change cell crystallized to various fractional levels. For these simulations, the microring is of diameter 19.05 µm, the gap between the ring and the coupling waveguide is 150 nm, and the rest of the waveguide and film dimensions are as given in Fig. 3. It is clear from the results of Fig. 7 that the readout signal (i.e., the optical transmission through the coupling waveguide) changes dramatically with the degree of crystallization. Indeed, as previously observed experimentally,² both the wavelength and the amplitude of the minimum in transmission change as the fraction of crystallized material changes (i.e., as the cell is set to various multilevel states).

Having demonstrated above some of the attributes and capabilities of our behavioral model, we now make a direct comparison of simulated results to experimental data. Thus, we wrote to a fully crystalline cell (on a straight Si waveguide device) with rectangular optical pulses of power 13.23 mW and a range of durations (specifically 10, 20, 50, 75, 100, 125, 150, 175, and 200 ns) while monitoring (using a continuous wave laser of low power) the waveguide’s optical transmission. The results are shown in Fig. 8, where a good match between experiment [Fig. 8(a)] and simulation [Fig. 8(b)] can be seen. In both cases (experiment and simulation), the shorter duration pulses do not induce any permanent change in transmission, meaning that little or no amorphization occurred. However, longer duration pulses (specifically 50 ns and above) did result in permanent changes in transmission, indicating successful amorphization. For pulse durations of 75 ns and above, the change in transmission saturated, indicating that, for this optical power and for this particular device, no significant increases in amorphized volume occurred in the cell for pulses longer than 75 ns.

In summary, we have developed a fast, computationally efficient behavioral model for integrated phase-change photonic devices which is capable of simulating their write, erase, and read processes. Thus, the model simulates the performance of such devices when used not only as binary and multilevel memories but also for arithmetic or logic processing and, importantly, as synaptic and neuronal mimics for neuromorphic computing applications. In this work, we have concentrated specifically on the use of the archetypal phase-change material GST, but the model is in no way limited simply to GST and could be used with any phase-change material so long as suitable material parameters were available (specifically, the real and imaginary parts of the refractive index, thermal conductivity and specific heat and density, along with activation energies, etc., used for the JMAK approach). Indeed, since the behavioral model runs in seconds as compared to several hours for equivalent FE models, it could provide a useful tool for the comparison of the likely performance of various phase-change alloys. Furthermore, it should prove extremely useful in the design and simulation of future integrated phase-change photonic systems, including neuromorphic photonic processors.

This work was supported via the European Union’s Horizon 2020 Research and Innovation Program (Grant No. 780848, Fun-COMP project) and the UK’s Engineering and Physical Sciences Research Council [the EPSRC CDT in Metamaterials (Grant No. EP/L015331/1), the WAFT Collaboration (Grant No. EP/M015173/1), and the Chalcogenide Advanced Manufacturing Partnership (Grant No. EP/M015130/1)]. Data relating to this paper available, on reasonable request, from the corresponding author.