NTMpy: An open source package for solving coupled parabolic differential equations in the framework of the three-temperature model
Lukas Alber, Valentino Scalera, Vivek Unikandanunni, Daniel Schick, Stefano Bonetti
NNTMpy: An open source package for solving coupled parabolic differential equations inthe framework of the three-temperature model
Lukas Alber a , Valentino Scalera b , Vivek Unikandanunni a , Daniel Schick c , Stefano Bonetti a,d a Department of Physics, Stockholm University, 106 91 Stockholm, Sweden b Department of Electrical Engineering and ICT, University of Naples Federico II, 80125 Naples, Italy. c Max Born Institute for Nonlinear Optics and Short Pulse Spectroscopy, Max-Born-Straße 2A, 12489 Berlin, Germany d Department of Molecular Sciences and Nanosystems, Ca’ Foscari University of Venice, 30172 Venezia-Mestre, Italy
Abstract
The NTMpy code package allows for simulating the one-dimensional thermal response of multilayer samples after opticalexcitation, as in a typical pump-probe experiment. Several Python routines are combined and optimized to solve coupledheat diffusion equations in one dimension, on arbitrary piecewise homogeneous material stacks, in the framework of theso-called three-temperature model. The energy source deposited in the material is modelled as a light pulse of arbitrarycross-section and temporal profile. A transfer matrix method enables the calculation of realistic light absorption inpresence 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 mathematicalbackground of the code, we lay out the workflow, and we validate the functionality of our package by comparing it tocommercial software, as well as to experimental transient reflectivity data recorded in a pump-probe experiment withfemtosecond light pulses.
Keywords:
Ultrafast condensed matter dynamics, Coupled parabolic differential equations, N-temperature model,Transfer Matrix Method
Program summary
Program title : NTMpy v.0.1.1
Program obtainable from : Github.com https://github.com/udcm-su/NTMpy
Licensing provisions : MIT license (MIT)
No. of lines in distributed program, including test data, etc. :1804
No. of bytes in distributed program, including test data, etc. :106212
Distribution format :(tar.gz)
Programming language : Python
Computer : PC, Workstation
Operating system : Running Python installation required. (Op-erating system independent)
RAM : Typically in the range of hundrets of megabytes
External routines : Python 3.5 or higher, numpy, matplotlib,bsplines, tqdm
Nature of problem : 1-dimensional coupled non linear partialdifferential equations; diffusion and relaxation dynamics for-multiple systems and multiple layers.
Solution method : Simulate the diffusion and relaxation dynam-ics of up to 3 coupled systems via an object oriented user inter-face. In order to approximate the solution and its derivatives inspace B-Spline interpolation is used. The solution is developedin time via the Explicit Euler method.
Unusual features : A routine to automatically select the idealtime step for stability of the algorithm is implemented. Rou- tines for output of raw data in order to post process and pre-made visualization routines are implemented.
Running time : Main influence factor is how many time stepsare required to integrate the equation in time. The time stepevaluation in the Explicit Euler loop depends on the input pa-rameters of the equation. Almost all tested problems can besolved under 30 seconds on a recent computer.
1. Introduction
The understanding of how heat is transported at thenanometer level and at ultrafast time scales in differentmaterials is an open question in modern condensed matterresearch. The increasing availability of commercial fem-tosecond laser systems makes it possible to study materialdynamics at sub-picosecond time scales, allowing for in-vestigating non-equilibrium energy transport. This alsoasks for numerical computations able to model the exper-imental evidence and to either validate or extract a set ofphysical parameters.A commonly used strategy to describe the highly non-equilibrium processes induced by ultrafast laser excitation,is the N -temperature model ( N = 2 or 3, typically) posedby Anisimov et al. [1]. This model is a set of coupledparabolic differential equations, which describe the tem-poral evolution of the energy of the electron, lattice, andspin systems as well as the transfer of energy between these Preprint submitted to Computer Physics Communications February 12, 2020 a r X i v : . [ c ond - m a t . o t h e r] F e b ystems. The N -temperature model is used to simulate abroad range of experiments in the ultrafast community,e.g. Ref. [2, 3, 4, 5], but an open source implementationof it is still missing.Here, we provide an open source, user-friendly Pythonpackage able to solve the N -temperature model in one-dimension, with arbitrary, piece-wise homogeneous layerswith different physical properties [6]. The time-step selec-tion to optimize the solving routines and guarantee theirconvergence, the calculation of the deposited energy byan arbitrarily shaped laser pulse using the transfer matrixmethod, and visualization routines, are all automatizedfor the early-user. However, the object oriented designof the package, also allows for easy customization for theadvanced users.After laying out the mathematical background on whichthe implementation of the solver is based on in Sec. 2,we introduce the workflow of the program by showing themost important commands and the output, users can ex-pect in Sec. 3. To validate our results we compare thecomputation of this open source package to commercialsoftware, and further more to real experimental data, inSec. 4.
2. Mathematical Methods and Background
When a material is illuminated by light, the energy ofthe electromagnetic radiation is partly reflected, partlytransmitted through the material, and partly absorbed byit. The absorbed electromagnetic energy essentially heatsthe material, rising its temperature. While this is a simpleproblem to solve for the case of a continuous light source,it becomes immediately more complex for the case of anultrashort laser pulses, i.e. with sub-picosecond pulse du-ration. In this case, the concept of temperature is in itselfan ill-defined one, because the different energy reservoirsin a material (the electronic E and lattice L systems, andalso the spin S system for a magnetically ordered sample)responds on different time scales and are not in equilib-rium among each other. Generally, only the electrons canreact fast enough to the ultrashort laser pulse, and theenergy is released from the electronic system to the otherheat 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 tosolve the heat equation in time, and then coupled to allowfor the energy to flow between them. This is addressed byconsidering three independent reservoirs, in order to allowfor the definition of a temperature and to solve the heatequation in time in each of the systems. Finally, the threesystems are coupled such that the energy can flow betweenthem. This idea is generalized to a N -temperature model,describing the heat exchange between the systems and theheat diffusion along the one-dimensional, multiple layeredmaterial. Each of the systems - electron, lattice, and spin- has its individual temperature T E,L,S ( x, t ) , where x de-notes the depth in the specimen and t is the time. The dynamics of every system are described by aparabolic partial differential equation, namely C E ( T E ) · ρ · ∂ t T E = ∂ x (cid:0) k E ( T E ) · ∂ x T E (cid:1) + G EL · ( T L − T E ) + G SE · ( T S − T E ) + S ( x, t ) ,C L ( T L ) · ρ · ∂ t T L = ∂ x (cid:0) k L ( T L ) · ∂ x T L (cid:1) + G EL · ( T E − T L ) + G LS · ( T S − T L ) ,C S ( T S ) · ρ · ∂ t T S = ∂ x (cid:0) k S ( T S ) · ∂ x T S (cid:1) + G SE · ( T E − T S ) + G LS · ( T L − T S ) . (1)Here C E,L,S ( T ) and k E,L,S ( T ) are the specific heat ca-pacity and thermal conductivity, which can be defined asconstants 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 is ( T i ) and k is ( T i ) , the dependence on the layer s will beomitted when not necessary in order to have a simpler no-tation.The term S ( x, t ) is responsible for the heat injection to theelectronic system and physically corresponds to a pulsedlaser source hitting the sample at the surface. Note thatin Eq. (1), we assume the coupling G , responsible for theheat exchange and relaxation between the systems, to belinear 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 subjectto research [9].The temperature profile is expressed as a linear com-bination of basis functions B m ( x ) for m = 1 , . . . , M with time depending coefficients for every subsystem, i ∈{ E, L, S } , c im ( t ) , i.e. T i ( x, t ) = M (cid:88) m =1 c im ( t ) B m ( x ) . (2)The solution of Eq. (1) is now reduced to a finite dimen-sional problem and consequently the diffusion equationcannot generally be solved exactly in the whole domain.The solution can be approximated by imposingEq. (1) to be satisfied on a given grid of points { x , x , x , . . . , x M − } , i.e. collocation points . Two ad-ditional points x and x M are added to the grid on theboundary of the domain to impose the boundary condi-tions.The temperatures and their derivatives have an exactanalytic expression at the grid points T i ( x j , t ) = M (cid:88) m =1 c im ( t ) B m ( x j ) ,∂T i ∂x ( x j , t ) = M (cid:88) m =1 c im ( t ) ∂B m ∂x ( x j ) ,∂ T i ∂x ( x j , t ) = M (cid:88) m =1 c im ( t ) ∂ B m ∂x ( x j ) . (3)2y introducing the M × M matrices D , D and D withgeneric elements { D } jm = B m ( x j ) , { D } jm = ∂B m ∂x ( x j ) , { D } jm = ∂ B m ∂x ( x j ) , (4)the temperatures and their derivatives can be obtained bymatrix products. By using the Leibnitz formula Eq. (1)can be the reformulated as follows ρC i D dd t c i = k i ( D c i ) · D c i + d k is d( D c i ) · ( D c i ) + S i ( x, t )+ (cid:88) k ∈{ E,L,S } k (cid:54) = i G ik ( D ( c k − c i )) , (5)for i = { E, L, S } , such that c i = { c i , c i , . . . , c iM } is thevector of coefficients relative to the i -th temperature, thedot ( · ) denotes the elementwise product of two vectors and ( a ) denotes the vector whose elements are the squares ofthe elements of the vector a .The time evolution of the coefficient vectors c i is evalu-ated using the explicit Euler formula, which reads D c i ( t + ∆ t ) = D c i ( t ) + ∆ tD d c i d t ( t ) , (6)where the time derivative is calculated as shown in Eq. (5)and S is present only for i = E .When an analytical formula for d k i / d T i is unavailable,this derivative can be computed numerically without intro-ducing significant errors since the conductivity is typicallya regular function. Equation (6) must be completed withthe initial- and boundary conditions at the left and rightend of the material under consideration.The boundary conditions in this software can be ofDirichlet type or of a modified Neumann type accordingto whether the temperature or the heat flux is assigned atthe boundaries. Hence we have the two following optionsfor the boundary conditionsDirichlet : (cid:40) ( D )[ 0 , :] c i = T BC ( x , t )( D )[ M, :] c i = T BC ( x M , t ) , orNeumann : (cid:40) k i ( T i ( x ))(( D c i ) D )[ 0 , :] c i = H BC ( x , t ) k i ( T i ( x M ))(( D c i ) D )[ M, :] c i = H BC ( x M , t ) , where the notation [ 0 , :] and [ M, :] indicates the firstand the last row of the matrix and T BC or H BC are theassigned boundary conditions.The basis functions chosen for the approximation of thesolution are B-Splines realized with the Cox-de Boor algo-rithm. These splines are continuous and have continuous derivatives up to a chosen order. A set of B-Spline basisfunctions, as used in the software, are depicted in Fig. 1,where a stack of two layers is under consideration.The use of C k -continuous functions is justified by thenature of the physical problem: the presence of a disconti-nuity in the temperature would cause an infinite heat flux,while a discontinuity in its first derivative would imply afinite heat flux into an infinitesimal control volume unlessthe conductivity is discontinuous.Discontinuities in the conductivity can be present whenthe specimen is made of two or more different materialsstacked together. In this case the heat flowing into theinterface between the s -th and the ( s + 1) -th layer, for the i -th temperature is H Interface = k i,s ( T i ) lim x → x − I ∂ x T i − k i,s +1 ( T i ) lim x → x + I ∂ x T i = 0 , (7)where the second subscript of k indicates the layer, x I isthe 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 theinterface a different set of B-Spline is considered for eachlayer (see Fig. 1), and the continuity of the temperatureand the condition Eq. (7), i.e. the conservation of the heatflow, are imposed at the interface.Notice that when the coefficients c im are determined thelimit and derivative appearing in Eq. (7) are analyticallycomputed. A crucial point for the precision, speed and stability ofthe code is the choice of an appropriate time step. Whenthis is not supplied by the user, the time step is auto-matically determined according to a criterion guarantee-ing stability on one side, but keeping the running time ofthe simulation as small as possible on the other side.For a linear N -temperature system with a single layerthe dynamical matrix is given by ρ diag (cid:26) k E C E , k L C L , k S C S (cid:27) ⊗ ( D − D ) + 1 ρ M ⊗ I N , (8) Depth of material S h a p e o f i ´ t h s p li n e f un c t i o n Figure 1: A set of order 5-continuous B-Splines for a two layer sys-tem, used to represent the solution in space, where each layer has alength of 1. M = − G EL + G SE C E G EL C E G SE C E G EL C L − G EL + G LS C L G LS C L G SE C S G LS C S − G LS + G SE C S (9)where I N is the N dimensional identity matrix and ⊗ de-notes the Kronecker product. In this system, to preservestability, the time step needs to satisfy the condition ∆ t < | λ | max , (10)where | λ | max is the eigenvalue of the dynamical matrixwith the largest absolute value. We note that, ∆ t dependson the input parameters C , k , G and on the desired spa-tial resolution, represented in D and D , but not on theheating 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 bycontemporaneous large k i ’s (thermal conductivities) andsmall C i ’s (specific heats), which produce large eigenvaluesin absolute value. To this purpose, in the worst scenariomatrix, the k i ’s are replaced by their maximum, evalu-ated on a set of values of T i and the C i ’s are replaced bytheir minimum evaluated on the same set. This yields thefollowing matrix ρ diag (cid:26) k Emax C Emin ρ , k
Lmax C Lmin ρ , k
Smax C Smin ρ (cid:27) ⊗ ( D − D )+ 1 ρ M max ⊗ I N , (11)where M max = − G EL + G SE C Emin G EL C Emin G SE C Emin G EL C Lmin − G EL + G LS C Lmin G LS C Lmin G SE C Smin G LS C Smin − G LS + G SE C Smin (12)As this linearization, from k i ( T ) to k maxi , depends onthe set of values chosen for T , typical values of temper-atures reached in the experiments on ultrafast dynam-ics are considered in the default case. That is, T ∈ [270,3000]K, but they can be changed by the user via sim.stability_lim([lowlimit,highlimit]) . Consid-ering the worst scenario allows to evaluate the time steponly once at the beginning of the simulation.In multilayer systems the above procedure is re-peated for each layer obtaining ∆ t , ∆ t , . . . , ∆ t L where L is the number of layers, the time step is ∆ t =min { ∆ t , ∆ t , . . . , ∆ t L } .
3. Background and Implementation
The code is implemented in Python with dependenceon the numpy and the bspline package for the numericalcomputation, on the matplotlib library for plotting of theresults and the progressbar package to monitor the elapsedtime. Installation of the NTMpy package will automati-cally check if the dependencies are fulfilled and downloadthe additional software if not.In order to make the user interface friendly, the code isobject oriented and it can be used either with commandline or in a script. The package contains three main classeswhich are source , simulation , and visual : While the source class is a collection of methods for the generation ofthe energy injection matrix, the simulation class handlesthe computation of the solution and builds the core of theprogram. After a solution has been found, it can be passedon to the visual classes allowing a fast and easy depictionof the results respectively. Even though in principle the source function S ( x, t ) , in-jecting heat in space and time, can be of generic type, themost common types of sources are already defined in thecode and one needs to specify only the main properties.That is, the user can choose between Lambert Beer’s lawor the transfer matrix method (TMM) to calculate theabsorption profile in space and independently from thatselect either a Gaussian, a repeated Gaussian or even acustom time profile to evaluate the shape of the heatingsource in time.For example a source with a Gaussian profile in timeand exponential decay in space, considering multiple re-flections, incident angle and polarization, i.e. TMM, isintroduced with the following lines of code. Note, that there is also a predefined way to calculatethe spacial absorption according to Lambert Beer´s law via s.spaceprofile = ’LB’ , which follows, except for the in-put of the incident angle and the polarization of the light,the same commands as above. Alternatively a custom pro-file in time can be given for the source, if arrays of data,4ere my_time and my_intensity , are provided to the pro-gram. Now considering the Lambert Beer decay law inspace, the commands to initialize a customized source are
After the initialization of a source, one can proceed toinitialize the simulation object, providing material specificparameters for every layer of the stack under considerationand then run the simulation. The constructor of the simu-lation class has a mandatory input which is the number N of temperatures under consideration and an optional inputwhich is the source s . N can be 1, 2 or 3. When the properties of s are not spec-ified as shown above, by default it is assumed that thereis no source term in Eq. (1), i.e. the fluence is 0.The properties of the media are assigned by adding se-quentially the layers of the materials with the method addLayer() . For instance where the inputs of the method are in order the lengthof the layer ( length ), the complex refractive index of thelayer ( n ), the list of the thermal conductivities of eachtemperature system ( [k1, k2] ), the list of the specificheats of each temperature system ( [C1, C2] ), the massdensity ( density ) and the exchange coupling ( G ).Thermal conductivities and specific heats can either benumbers or lambda functions when they vary with tem-perature. In case N = 3 , G can be a list containing thecoupling constants G , G and G in this order.Before running the simulation, a final time must be givenby the user through the command The commands described so far are sufficient to run thesimulation. If no further input is given some default valuesare selected according to the most common conditions ormodels employed for the experiments.By default the time step is assigned automatically throughthe procedure illustrated in Sec. 2.1. However the user candefine a different time step by sim . time_step = time_step
The usual initial condition for all the temperatures T i bydefault is T i = 300 K. The user can modify the initial con-dition through the command line where i varies between 1 and N , indicating, that differentinitial conditions can be set for each subsystem. Further-more the initial temperature T0 can be either a numberor a lambda function describing the space profile of thetemperature T ( x, .By default the boundary conditions are of modified Neu-mann type with no heat flux through the boundaries, cor-responding to an insulated system. The boundary condi-tion type and value can be modified through the followingmethods where again i is the index of the temperature system sub-ject to the change, side can be either ’left’ or ’right’ depending on which side the change is applied to, Type iseither ’dirichlet’ or ’neumann’ respectively, and BC isa constant or a lambda function providing the value of theboundary condition as the time varies.Once the input of the data is complete the simulation isexecuted by the command It yields the numpy.array containing the 2+N-dimensionalarrays T of the temperature of the N subsystems on thespace grid x at the times t . Once the simulation has been executed the results, thatis the dynamics of all respective temperatures in space andtime, are provided to the user in array form, as describedabove. However, in order to make the visualization both,easier and quicker for the the user, a separate class withrelevant plotting methods has been defined. This givesthe user the freedom to output the raw data and do post-processing on their own but also to use the visualizationclass in order to depict some properties immediately.Among others, there are methods to visualize the heat-ing source S ( x, t ) , a contour plot of the temperature inspace and time T E,L,S ( x, t ) , but also some simple readymade post-process routines like the plot of the averagetemperature of a layer over time, see Fig. 2. Moreover thepackage provides a method to play an animation of how allthe temperature systems evolve in time which can be use-ful either to visually detect the impact of the fundamentalmechanisms in the case considered, or for pedagogical use.Different routines can be called by following the commands Time (ps) T e m p e r a t u r e ( K ) System 1System 2System 30.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0Distance from surface (nm)02468 T i m e ( p s ) T e m p e r a t u r e ( K ) a)b) Figure 2: Example of a 3-temperature simulation and the output oftwo visualization methods: a) v.average() and b) v.contour(’1’) described in the text. Here a
20 nm ferromagnetic nickel thin film,heated by a femtosecond laser source is under consideration. Thephysical parameters are taken from Ref. [2]. where the output of the first two visualization routines arearrays, which allow users, to look and process the raw datathemselves. The only visualization methods that requireinput arguments are v.contour(’N’) , where N can be 1,2or 3, corresponding to electron-, the lattice- or the spin-system, and animation(speed,save) , where the speed ofthe animation can be adjusted and the user can decidewhether they want to save the animation, with save = 1 or not, with save = 0 .This makes it easy to quickly check results and analyzeproperties of the dynamics of the simulation.
4. Package Validation
In order to validate the results obtained from theNTMpy package, the focus is set on the two main parts ofthe software. First, the evaluation of the local absorption with respect to the incident laser pulse, and second, theperformance of the solver itself. Looking at Eq. (1), we areinterested of how accurately we are computing the sourceterm S ( x, t ) responsible for heating and how well our ob-tained solution T ( x, t ) matches with commercial software,when experimentally relevant cases are under investiga-tion.For the evaluation of the local absorption of S ( x, t ) , weare considering a [Pt 3nm|Co 15nm|Cr 5nm|MgO]- mate-rial stack and compute the local absorption per unit in-cident power, with respect to the distance from the sur-face, with the implemented transfer matrix method. Theobtained profile is compared to the result, obtained fromCOMSOL Multiphysics™, Ref. [10], simulations, wherethe same light and material parameters are considered. InFig. 3 a), we see, that both, the commercial and our opensource software show an overlapping result (mean relativediscrepancy δ ( x ) = 2 . ).For the evaluation of computation time and accuracyof 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 pa-rameters for both software and obtain output that is verymuch in agreement, with respect to space Fig. 3 b) andtime c) dimension. That is the relative error of the temper-ature enhancement T NTMpy − T Matlab T Matlab − ≤ . Also the timeneeded until a solution is obtained from our free packageis comparable to Matlab’s performance. Even though the software is able to solve generic couplednon linear 1-D diffusion equations in the form of Eq. (1),the focus is laid on heat transfer between systems andtransport along the material, within the framework of thetwo/three temperature model.If the physical parameters of all the material layers un-der consideration are known, the software can be used todepict the temperature dynamics in space and time of theselected configuration. This provides an opportunity to beable to detect local phenomena in a material, as opposedto averaged effects, measured in most experiments. How-ever reliable data on certain parameters is sparse and it isa current interest of research to find those properties.In order to demonstrate the use of this software, we mea-sured the change in reflectivity on a
10 nm platinum thinfilm on a silicon substrate and used the NTMpy software toperform a simulation of the temperature dynamics of thesame system. The experimental data was retrieved from apump-probe transient reflectivity measurement. Accord-ing to Refs. [7, 12], a linear relationship between the mea-sured change in reflectivity and the change in temperaturecan be assumed, i.e. ∆ RR ∝ ∆ TT . This makes it possible tocompare the experiment to the simulation.In the simulation, we used the same nominal laser sourceparameters as in the experiment. That is a
400 nm , p-polarized laser pulse with a FWHM of
100 fs and a fluence6
Distance from surface (nm) A b s o r p t i o n p e r i n c i d e n t p o w e r ( n m ) NTMpyInterpolated COMSOL data
Distance from surface (nm) T e m p e r a t u r e ( K ) NTMpyMatlab
Time (ps) T e m p e r a t u r e ( K ) NTMpyMatlab a) b) c)
Figure 3: (a) Local absorption per incident laser power computed with NTMpy’s TMM-module and with the commercial software COMSOLMultiphysics™, Ref. [10]. For both simulations a three-layer material with identical parameters is considered. Line cuts for two coupledtemperature systems T E,L ( x, t ) at (b) a fixed time t = t and at (c) a fixed point in space x = 0 . NTMpy’s solutions for both systems(electronic and lattice), computed for a two-layer stack, are compared to Matlab’s built-in partial differential equation pdepe() solver, Ref.[11].Figure 4: Local power absorption for the platinum/silicon systemunder consideration, calculated via the transfer matrix method im-plemented in NTMpy. Here a p-polarized laser pulse with a wave-length of
400 nm , hitting the sample at a ° angle with a fluence of / cm and a FWHM of . is considered. of / cm of which 53% get absorbed, which is hittingthe sample in a 45° angle. The local absorption profile isshown in Fig. 4. For the platinum film and the siliconsubstrate the parameters from Table 1 are used.The simulation is set up as follows: sim = simulation (2 ,s)sim . addLayer ( length_Pt ,n_Pt ,[ k_el , k_lat ],[C_el , C_lat ],rho ,[ G_Pt ])sim . addSubstrate ("Si")sim . final_time = 7*u.ps[x,t,T] = sim . run () Time (ps) N o r m a li z e d t e m p e r a t u r e / r e f l e c t i v i t y Simulation with substrateSimulation without substrateExperiment
Figure 5: Experimental data of a reflectivity measurement (symbols)compared to the simulated temperature dynamics of a platinum thinfilm on a silicon substrate (solid line) and free-standing (dashed line).A linear relation between the change in temperature and change inreflectivity is assumed, and they are normalized to the maximumchange. The simulation parameters for the materials are taken fromTable 1.
To compare the simulated data with the experimentalcase, we computed the exponentially weighted tempera-ture data T E,L ( x, t ) with respect to the depth of the ma-terial. This mimics the effect of limited optical penetrationdepth of the probe laser. In Fig. 5, we show a comparisonbetween a simulation with and without the silicon sub-strate.The results, depicted in Fig. 5 show that the multilayersolution obtained from the NTMpy software is able to re-produce the experimental data to a very good degree. Fur-7arameter Symbol/ Units Platinum SiliconThickness l (nm) 10 100000Refractive index (at 400 nm) n k e (Wm − K − ) 72 [15] 130 [16]Lattice heat conductivity k l (Wm − K − ) 72 [15] k l ( T ) [17]Electron heat capacity C e ( T )(Jkg − K − ) /ρ P t T e [15] 150/ ρ Si · T e [17]Lattice heat capacity C l (Jkg − K − ) 2.78E6/ ρ P t [15] 1.6E6/ ρ Si [17]Electron phonon coupling G (Wm − K − ) Table 1: Material parameters used for the TTM simulation in Fig. 5. thermore, one can clearly see that disregarding the sub-strate in the simulation, which is the simplest way to solvethe two-temperature model, leads to a considerably dif-ferent results. In this case, a four times larger electron-lattice coupling constant from G = 2 . · (Wm − K − ) to G ≈ · (Wm − K − ) would be needed to reestab-lish agreement with the experiment.
5. Conclusion
We implemented NTMpy, an open source Python basedsoftware package for solving coupled parabolic differentialequations in one dimension. Along with the mathematicalbackground of the algorithm, we introduced the structureand the work flow of the program. The object orientedway in which the program is designed gives the user thefreedom to run tailor-made simulations, without havingto worry about coding details, which are automatized orready-made. This helps to focus more on the output andanalysis of the fundamental dynamics. Together with thevisualization class, the software can not only be used by re-searchers to investigate relaxation and diffusion dynamics,but also for educational purpose. Finally, a comparisonof the output of NTMpy with the one of both commercialsoftware and also of new experimental data, demonstratedthe numerical reliability of the software, and its ability toproduce physically realistic results.
Acknowledgements
L.A. acknowledges support from the Swedish ResearchCouncil (VR), Grant: 2018-04611. V.U. and S.B. acknowl-edge support from the European Research Council, Start-ing Grant 715452 “MAGNETIC-SPEED-LIMIT”.
References [1] B. L. K. S. I. Anisimov, T. L. Perel’man, Electron emission frommetal surfaces exposed to ultrashort laser pulses, Sov. Phys.JETP 39 (1974) 375–377 ((February 1974).[2] E. Beaurepaire, J.-C. Merle, A. Daunois, J.-Y. Bigot,Ultrafast spin dynamics in ferromagnetic nickel,Phys. Rev. Lett. 76 (1996) 4250–4253 (May 1996). doi:10.1103/PhysRevLett.76.4250 .URL https://link.aps.org/doi/10.1103/PhysRevLett.76.4250 [3] A. A. Maznev, J. A. Johnson, K. A. Nelson, Non-equilibriumtransient thermal grating relaxation in metal, Journal of Ap-plied Physics 109 (7) (2011) 073517 (2011). arXiv:https://doi.org/10.1063/1.3569731 , doi:10.1063/1.3569731 .URL https://doi.org/10.1063/1.3569731 [4] B. Koopmans, G. Malinowski, F. Dalla Longa, D. Steiauf,M. Fähnle, T. Roth, M. Cinchetti, M. Aeschlimann, Explainingthe paradoxical diversity of ultrafast laser-induced demagneti-zation, Nature materials 9 (3) (2010) 259 (2010).[5] D. Schick, A. Bojahr, M. Herzog, R. Shayduk, C. Von Ko-rff Schmising, M. Bargheer, UDKM1DSIM - A simulationtoolkit for 1D ultrafast dynamics in condensed matter, Com-puter Physics Communications 185 (2) (2014) 651–660 (2014). doi:10.1016/j.cpc.2013.10.009 .[6] L. Alber, V. Scalera, V. Unikandanunni, S. Bonetti, udcm-su/ntmpy: Ntmpy-package1.1.1 (2019). doi:10.5281/ZENODO.3553742 .URL https://zenodo.org/record/3553742 [7] J. Hohlfeld, S.-S. Wellershoff, J. Güdde, U. Conrad, V. Jähnke,E. Matthias, Electron and lattice dynamics following opticalexcitation of metals, Chemical Physics 251 (2000) 237–258 (012000). doi:10.1016/S0301-0104(99)00330-4 .[8] P. M. Norris, A. P. Caffrey, R. J. Stevens, J. M. Klopf, J. T.McLeskey, A. N. Smith, Femtosecond pump–probe nondestruc-tive examination of materials (invited), Review of Scientific In-struments 74 (1) (2003) 400–406 (2003). arXiv:https://doi.org/10.1063/1.1517187 , doi:10.1063/1.1517187 .URL https://doi.org/10.1063/1.1517187 [9] Z. Lin, L. V. Zhigilei, V. Celli, Electron-phonon coupling andelectron heat capacity of metals under conditions of strongelectron-phonon nonequilibrium, Phys. Rev. B 77 (2008)075133 (Feb 2008). doi:10.1103/PhysRevB.77.075133 .URL https://link.aps.org/doi/10.1103/PhysRevB.77.075133 [10] COMSOL Multiphysics™, , COMSOL AB,v.5.4.[11] MATLAB, , The MathWorks Inc., v. R2018b(9.5.0.944444).[12] A. P. Caffrey, P. E. Hopkins, J. M. Klopf, P. M. Norris,Thin film non-noble transition metal thermophysical properties,Microscale Thermophysical Engineering 9 (4) (2005) 365–377(2005). arXiv:https://doi.org/10.1080/10893950500357970 , doi:10.1080/10893950500357970 .URL https://doi.org/10.1080/10893950500357970 [13] A. D. Rakić, A. B. Djurišić, J. M. Elazar, M. L. Majewski,Optical properties of metallic films for vertical-cavity optoelec-tronic devices, Appl. Opt. 37 (22) (1998) 5271–5283 (Aug 1998). doi:10.1364/AO.37.005271 .URL http://ao.osa.org/abstract.cfm?URI=ao-37-22-5271 [14] D. E. Aspnes, A. A. Studna, Dielectric functions and opticalparameters of si, ge, gap, gaas, gasb, inp, inas, and insb from1.5 to 6.0 ev, Phys. Rev. B 27 (1983) 985–1009 (Jan 1983). doi:10.1103/PhysRevB.27.985 .URL https://link.aps.org/doi/10.1103/PhysRevB.27.985 [15] D. E. Gray, I. E. Dayton, American institute of physics hand-book, American Journal of Physics 32 (5) (1964) 389–389(1964). arXiv:https://doi.org/10.1119/1.1970399 , doi:10. .URL https://doi.org/10.1119/1.1970399 [16] C. J. Glassbrenner, G. A. Slack, Thermal conductivity of siliconand germanium from 3°k to the melting point, Phys. Rev. 134(1964) A1058–A1069, taken at room temperature (May 1964). doi:10.1103/PhysRev.134.A1058 .URL https://link.aps.org/doi/10.1103/PhysRev.134.A1058 [17] I. Lazanu, S. Lazanu, Contribution of the electron–phonon in-teraction to lindhard energy partition at low energy in ge andsi detectors for astroparticle physics applications, AstroparticlePhysics 75 (2016) 44 – 54 (2016). doi:https://doi.org/10.1016/j.astropartphys.2015.09.007 .URL [18] J. Hohlfeld, S.-S. Wellershoff, J. Güdde, U. Conrad, V. Jähnke,E. Matthias, Electron and lattice dynamics following opti-cal excitation of metals, Chemical Physics 251 (1) (2000) 237– 258 (2000). doi:https://doi.org/10.1016/S0301-0104(99)00330-4 .URL [19] A. Chettah, H. Kucal, Z. Wang, M. Kac, A. Meftah, M. Toule-monde, Behavior of crystalline silicon under huge electronic ex-citations: A transient thermal spike description, Nuclear In-struments and Methods in Physics Research Section B: BeamInteractions with Materials and Atoms 267 (16) (2009) 2719– 2724, proceedings of the 23rd International Conference onAtomic Collisions in Solids (2009). doi:https://doi.org/10.1016/j.nimb.2009.05.063 .URL