NTMpy: An open source package for solving coupled parabolic differential equations in the framework of the three-temperature model

The NTMpy code package allows for simulating the one-dimensional thermal response of multilayer samples after optical excitation, as in a typical pump-probe experiment. Several Python routines are combined and optimized to solve coupled heat diffusion equations in one dimension, on arbitrary piecewise homogeneous material stacks, in the framework of the so-called three-temperature model. The energy source deposited in the material is modelled as a light pulse of arbitrary cross-section and temporal profile. A transfer matrix method enables the calculation of realistic light absorption in presence of scattering interfaces as in multilayer samples. The open source code is fully object-oriented to enable a user-friendly and intuitive interface for adjusting the physically relevant input parameters. Here, we describe the mathematical background of the code, we lay out the workflow, and we validate the functionality of our package by comparing it to commercial software, as well as to experimental transient reflectivity data recorded in a pump-probe experiment with femtosecond light pulses.


Introduction
The understanding of how heat is transported at the nanometer level and at ultrafast time scales in different materials is an open question in modern condensed matter research. The increasing availability of commercial femtosecond laser systems makes it possible to study material dynamics at sub-picosecond time scales, allowing for investigating non-equilibrium energy transport. This also asks for numerical computations able to model the experimental evidence and to either validate or extract a set of physical parameters.
A commonly used strategy to describe the highly nonequilibrium processes induced by ultrafast laser excitation, is the N -temperature model (N = 2 or 3, typically) posed by Anisimov et al. [1]. This model is a set of coupled parabolic differential equations, which describe the temporal evolution of the energy of the electron, lattice, and spin systems as well as the transfer of energy between these systems. The N -temperature model is used to simulate a broad range of experiments in the ultrafast community, e.g. Ref. [2,3,4,5], but an open source implementation of it is still missing.
Here, we provide an open source, user-friendly Python package able to solve the N -temperature model in onedimension, with arbitrary, piece-wise homogeneous layers with different physical properties [6]. The time-step selection to optimize the solving routines and guarantee their convergence, the calculation of the deposited energy by an arbitrarily shaped laser pulse using the transfer matrix method, and visualization routines, are all automatized for the early-user. However, the object oriented design of the package, also allows for easy customization for the advanced users.
After laying out the mathematical background on which the implementation of the solver is based on in Sec. 2, we introduce the workflow of the program by showing the most important commands and the output, users can expect in Sec. 3. To validate our results we compare the computation of this open source package to commercial software, and further more to real experimental data, in Sec. 4.

Mathematical Methods and Background
When a material is illuminated by light, the energy of the electromagnetic radiation is partly reflected, partly transmitted through the material, and partly absorbed by it. The absorbed electromagnetic energy essentially heats the material, rising its temperature. While this is a simple problem to solve for the case of a continuous light source, it becomes immediately more complex for the case of an ultrashort laser pulses, i.e. with sub-picosecond pulse duration. In this case, the concept of temperature is in itself an ill-defined one, because the different energy reservoirs in a material (the electronic E and lattice L systems, and also the spin S system for a magnetically ordered sample) responds on different time scales and are not in equilibrium among each other. Generally, only the electrons can react fast enough to the ultrashort laser pulse, and the energy is released from the electronic system to the other heat reservoirs only at a later time.
To treat this problem, three reservoirs are considered, in order to allow for the definition of a temperature and to solve the heat equation in time, and then coupled to allow for the energy to flow between them. This is addressed by considering three independent reservoirs, in order to allow for the definition of a temperature and to solve the heat equation in time in each of the systems. Finally, the three systems are coupled such that the energy can flow between them. This idea is generalized to a N -temperature model, describing the heat exchange between the systems and the heat diffusion along the one-dimensional, multiple layered material. Each of the systems -electron, lattice, and spin -has its individual temperature T E,L,S (x, t), where x denotes the depth in the specimen and t is the time.
The dynamics of every system are described by a parabolic partial differential equation, namely Here C E,L,S (T ) and k E,L,S (T ) are the specific heat capacity and thermal conductivity, which can be defined as constants or as functions of the respective temperature T E,L,S . Since we consider a stack of multiple layers C i and k i also depend on the respective layer, hence they are C i s (T i ) and k i s (T i ), the dependence on the layer s will be omitted when not necessary in order to have a simpler notation. The term S(x, t) is responsible for the heat injection to the electronic system and physically corresponds to a pulsed laser source hitting the sample at the surface. Note that in Eq. (1), we assume the coupling G, responsible for the heat exchange and relaxation between the systems, to be linear as in the formulation from Anisimov et al., Ref. [7], [8] and other groups. However, it should be mentioned, that considering a non-constant G(T ) is current subject to research [9].
The temperature profile is expressed as a linear combination of basis functions B m (x) for m = 1, . . . , M with time depending coefficients for every subsystem, i ∈ {E, L, S}, c i m (t), i.e.
The solution of Eq. (1) is now reduced to a finite dimensional problem and consequently the diffusion equation cannot generally be solved exactly in the whole domain. The solution can be approximated by imposing Eq. (1) to be satisfied on a given grid of points {x 1 , x 2 , x 3 , . . . , x M −1 }, i.e. collocation points. Two additional points x 0 and x M are added to the grid on the boundary of the domain to impose the boundary conditions.
The temperatures and their derivatives have an exact analytic expression at the grid points By introducing the M × M matrices D 0 , D 1 and D 2 with generic elements the temperatures and their derivatives can be obtained by matrix products. By using the Leibnitz formula Eq. (1) can be the reformulated as follows . . , c i M } is the vector of coefficients relative to the i-th temperature, the dot (·) denotes the elementwise product of two vectors and (a) 2 denotes the vector whose elements are the squares of the elements of the vector a.
The time evolution of the coefficient vectors c i is evaluated using the explicit Euler formula, which reads where the time derivative is calculated as shown in Eq. (5) and S is present only for i = E.
When an analytical formula for dk i /dT i is unavailable, this derivative can be computed numerically without introducing significant errors since the conductivity is typically a regular function. Equation (6) must be completed with the initial-and boundary conditions at the left and right end of the material under consideration.
The boundary conditions in this software can be of Dirichlet type or of a modified Neumann type according to whether the temperature or the heat flux is assigned at the boundaries. Hence we have the two following options for the boundary conditions Dirichlet : where the notation  Fig. 1, where a stack of two layers is under consideration.
The use of C k -continuous functions is justified by the nature of the physical problem: the presence of a discontinuity in the temperature would cause an infinite heat flux, while a discontinuity in its first derivative would imply a finite heat flux into an infinitesimal control volume unless the conductivity is discontinuous.
Discontinuities in the conductivity can be present when the specimen is made of two or more different materials stacked together. In this case the heat flowing into the interface between the s-th and the (s + 1)-th layer, for the i-th temperature is where the second subscript of k indicates the layer, x I is the position of the interface and the superscript + and − indicates the limit from right and the left respectively.
In order to correctly represent the temperatures at the interface a different set of B-Spline is considered for each layer (see Fig. 1), and the continuity of the temperature and the condition Eq. (7), i.e. the conservation of the heat flow, are imposed at the interface.
Notice that when the coefficients c im are determined the limit and derivative appearing in Eq. (7) are analytically computed.

Evaluation of the Time Step
A crucial point for the precision, speed and stability of the code is the choice of an appropriate time step. When this is not supplied by the user, the time step is automatically determined according to a criterion guaranteeing stability on one side, but keeping the running time of the simulation as small as possible on the other side.
For a linear N -temperature system with a single layer the dynamical matrix is given by    9) where I N is the N dimensional identity matrix and ⊗ denotes the Kronecker product. In this system, to preserve stability, the time step needs to satisfy the condition where |λ| max is the eigenvalue of the dynamical matrix with the largest absolute value. We note that, ∆t depends on the input parameters C, k, G and on the desired spatial resolution, represented in D 0 and D 2 , but not on the heating source term S(x, t).
In the nonlinear context we can still use condition (10) with the eigenvalues of an adapted version of matrix (8). We consider the worst scenario for stability represented by contemporaneous large k i 's (thermal conductivities) and small C i 's (specific heats), which produce large eigenvalues in absolute value. To this purpose, in the worst scenario matrix, the k i 's are replaced by their maximum, evaluated on a set of values of T i and the C i 's are replaced by their minimum evaluated on the same set. This yields the following matrix In multilayer systems the above procedure is repeated for each layer obtaining ∆t 1 , ∆t 2 , . . . , ∆t L where L is the number of layers, the time step is ∆t = min{∆t 1 , ∆t 2 , . . . , ∆t L }.

Background and Implementation
The code is implemented in Python with dependence on the numpy and the bspline package for the numerical computation, on the matplotlib library for plotting of the results and the progressbar package to monitor the elapsed time. Installation of the NTMpy package will automatically check if the dependencies are fulfilled and download the additional software if not.
In order to make the user interface friendly, the code is object oriented and it can be used either with command line or in a script. The package contains three main classes which are source, simulation, and visual: While the source class is a collection of methods for the generation of the energy injection matrix, the simulation class handles the computation of the solution and builds the core of the program. After a solution has been found, it can be passed on to the visual classes allowing a fast and easy depiction of the results respectively.

Data Input
Even though in principle the source function S(x, t), injecting heat in space and time, can be of generic type, the most common types of sources are already defined in the code and one needs to specify only the main properties. That is, the user can choose between Lambert Beer's law or the transfer matrix method (TMM) to calculate the absorption profile in space and independently from that select either a Gaussian, a repeated Gaussian or even a custom time profile to evaluate the shape of the heating source in time.
For example a source with a Gaussian profile in time and exponential decay in space, considering multiple reflections, incident angle and polarization, i.e. TMM, is introduced with the following lines of code. Note, that there is also a predefined way to calculate the spacial absorption according to Lambert Beer s law via s.spaceprofile = 'LB', which follows, except for the input of the incident angle and the polarization of the light, the same commands as above. Alternatively a custom profile in time can be given for the source, if arrays of data, here my_time and my_intensity, are provided to the program. Now considering the Lambert Beer decay law in space, the commands to initialize a customized source are # Set source type s . spaceprofile = ' LB ' s . timeprofile = ' Custom ' # Set the value for the source s . loadData = [ my_time , my_intensity ] After the initialization of a source, one can proceed to initialize the simulation object, providing material specific parameters for every layer of the stack under consideration and then run the simulation. The constructor of the simulation class has a mandatory input which is the number N of temperatures under consideration and an optional input which is the source s. where the inputs of the method are in order the length of the layer (length), the complex refractive index of the layer (n), the list of the thermal conductivities of each temperature system ([k1, k2]), the list of the specific heats of each temperature system ([C1, C2]), the mass density (density) and the exchange coupling (G). Thermal conductivities and specific heats can either be numbers or lambda functions when they vary with temperature. In case N = 3, G can be a list containing the coupling constants G 12 , G 23 and G 31 in this order.
Before running the simulation, a final time must be given by the user through the command The usual initial condition for all the temperatures T i by default is T i = 300K. The user can modify the initial condition through the command line # Modify initial condition of i -th temp . ( in K ) sim . changeInit ( i , T0 ) where i varies between 1 and N, indicating, that different initial conditions can be set for each subsystem. Furthermore the initial temperature T0 can be either a number or a lambda function describing the space profile of the temperature T (x, 0).
By default the boundary conditions are of modified Neumann type with no heat flux through the boundaries, corresponding to an insulated system. The boundary condition type and value can be modified through the following methods # Change BC type for the i -th temp sim . changeBC_Type ( i , side , Type ) # Change BC value for the i -th temp sim . c hangeB C_Valu e ( i , side , BC ) where again i is the index of the temperature system subject to the change, side can be either 'left' or 'right' depending on which side the change is applied to, Type is either 'dirichlet' or 'neumann' respectively, and BC is a constant or a lambda function providing the value of the boundary condition as the time varies.
Once the input of the data is complete the simulation is executed by the command # Compute temperature map [x , t , T ] = sim . run () It yields the numpy.array containing the 2+N-dimensional arrays T of the temperature of the N subsystems on the space grid x at the times t.

Plot Results
Once the simulation has been executed the results, that is the dynamics of all respective temperatures in space and time, are provided to the user in array form, as described above. However, in order to make the visualization both, easier and quicker for the the user, a separate class with relevant plotting methods has been defined. This gives the user the freedom to output the raw data and do postprocessing on their own but also to use the visualization class in order to depict some properties immediately.
Among others, there are methods to visualize the heating source S(x, t), a contour plot of the temperature in space and time T E,L,S (x, t), but also some simple ready made post-process routines like the plot of the average temperature of a layer over time, see Fig. 2. Moreover the package provides a method to play an animation of how all the temperature systems evolve in time which can be useful either to visually detect the impact of the fundamental mechanisms in the case considered, or for pedagogical use. Different routines can be called by following the commands where the output of the first two visualization routines are arrays, which allow users, to look and process the raw data themselves. The only visualization methods that require input arguments are v.contour('N'), where N can be 1,2 or 3, corresponding to electron-, the lattice-or the spinsystem, and animation(speed,save), where the speed of the animation can be adjusted and the user can decide whether they want to save the animation, with save = 1 or not, with save = 0. This makes it easy to quickly check results and analyze properties of the dynamics of the simulation.

Comparison With Other Simulation Tools
In order to validate the results obtained from the NTMpy package, the focus is set on the two main parts of the software. First, the evaluation of the local absorption with respect to the incident laser pulse, and second, the performance of the solver itself. Looking at Eq. (1), we are interested of how accurately we are computing the source term S(x, t) responsible for heating and how well our obtained solution T (x, t) matches with commercial software, when experimentally relevant cases are under investigation.
For the evaluation of the local absorption of S(x, t), we are considering a [Pt 3nm|Co 15nm|Cr 5nm|MgO]-material stack and compute the local absorption per unit incident power, with respect to the distance from the surface, with the implemented transfer matrix method. The obtained profile is compared to the result, obtained from COMSOL Multiphysics™, Ref. [10], simulations, where the same light and material parameters are considered. In Fig. 3 a), we see, that both, the commercial and our open source software show an overlapping result (mean relative discrepancy δ(x) = 2.2%).
For the evaluation of computation time and accuracy of the obtained solution, we compare T (x, t) to Matlab s, Ref. [11], built in partial differential equation, pdepe(),solver. Again we consider the same light and material parameters for both software and obtain output that is very much in agreement, with respect to space Fig. 3 b) and time c) dimension. That is the relative error of the temperature enhancement Also the time needed until a solution is obtained from our free package is comparable to Matlab's performance.

Comparison With Experiment
Even though the software is able to solve generic coupled non linear 1-D diffusion equations in the form of Eq. (1), the focus is laid on heat transfer between systems and transport along the material, within the framework of the two/three temperature model.
If the physical parameters of all the material layers under consideration are known, the software can be used to depict the temperature dynamics in space and time of the selected configuration. This provides an opportunity to be able to detect local phenomena in a material, as opposed to averaged effects, measured in most experiments. However reliable data on certain parameters is sparse and it is a current interest of research to find those properties.
In order to demonstrate the use of this software, we measured the change in reflectivity on a 10 nm platinum thin film on a silicon substrate and used the NTMpy software to perform a simulation of the temperature dynamics of the same system. The experimental data was retrieved from a pump-probe transient reflectivity measurement. According to Refs. [7,12], a linear relationship between the measured change in reflectivity and the change in temperature can be assumed, i.e. ∆R R ∝ ∆T T . This makes it possible to compare the experiment to the simulation.
In the simulation, we used the same nominal laser source parameters as in the experiment. That is a 400 nm, ppolarized laser pulse with a FWHM of 100 fs and a fluence  of 6 mJ/cm 2 of which 53% get absorbed, which is hitting the sample in a 45°angle. The local absorption profile is shown in Fig. 4. For the platinum film and the silicon substrate the parameters from Table 1  To compare the simulated data with the experimental case, we computed the exponentially weighted temperature data T E,L (x, t) with respect to the depth of the material. This mimics the effect of limited optical penetration depth of the probe laser. In Fig. 5, we show a comparison between a simulation with and without the silicon substrate.
thermore, one can clearly see that disregarding the substrate in the simulation, which is the simplest way to solve the two-temperature model, leads to a considerably different results. In this case, a four times larger electronlattice coupling constant from G = 2.5 · 10 17 (Wm −3 K −1 ) to G ≈ 11 · 10 17 (Wm −3 K −1 ) would be needed to reestablish agreement with the experiment.

Conclusion
We implemented NTMpy, an open source Python based software package for solving coupled parabolic differential equations in one dimension. Along with the mathematical background of the algorithm, we introduced the structure and the work flow of the program. The object oriented way in which the program is designed gives the user the freedom to run tailor-made simulations, without having to worry about coding details, which are automatized or ready-made. This helps to focus more on the output and analysis of the fundamental dynamics. Together with the visualization class, the software can not only be used by researchers to investigate relaxation and diffusion dynamics, but also for educational purpose. Finally, a comparison of the output of NTMpy with the one of both commercial software and also of new experimental data, demonstrated the numerical reliability of the software, and its ability to produce physically realistic results.