Multilayer Microscopic Three Temperature Model Simulations (ML_M3TM_Sim)

Welcome! This code solves the M3TM equations in a magnetic material in a multilayer of an arbitrary amount of cap and substrate layers. This code is not yet an official python package, so to use it you must download the code and save it locally. This is a quick guide to set up a simulation, e.g. how to write a Scriptfile, specifying and shortly explaining the input parameters. The last version control has been done for Python 3.11.

Necessary packages

The following packages need to be installed to run the code:

• numpy
• scipy
• os
• time
• datetime
• matplotlib

Creating the Script file

Create a script within the code/Scripts folder. For here, we shall call our file tutorial_input.py. All the contents of this file are listed below.

Source code imports

How to import the source code

These are relative imports of all the source files needed for the simulation:

from ..Source.mats import SimMaterials
from ..Source.sample import SimSample
from ..Source.pulse import SimPulse
from ..Source.mainsim import SimDynamics

Create materials

How to create a material

Materials are defined merely by their parameters. Here is a list of all the available parameters, where [optional_1] denotes that these parameters need to be defined for materials with an electronic subsystem and [optional_2] parameters need to be introduced for materials with a spin subsystem:

• name (string). Name of the material
• tdeb (float). Debye temperature of the material
• cp_max (float). Maximal phononic heat capacity in W/m**3/K. Temperature dependence is computed with Einstein model
• kappap (float). Phononic heat diffusion constant in W/m/K
• kappae [optional_1] (float). Electronic heat diffusion constant in W/m/K
• ce_gamma [optional_1] (float). Sommerfeld constant of electronic heat capacity in J/m**3/K
• gep [optional_1] (float). Electron-phonon coupling constant in W/m**3/K
• spin [optional_2] (float). Effective spin of the material
• tc [optional_2] (float). Curie temperature of the material
• muat [optional_2] (float). Atomic magnetic moment in unit of \mu_Bohr
• asf [optional_2] (float). Electron-phonon-scattering induced spin flip probability of the material
• vat [optional_2] (float). Magnetic atomic volume in m^3.

For an insulating material, you only need to define the parameters for a phononic system:

my_insulator = SimMaterials(name='Gummydummy', tdeb=400, cp_max=3e6, kappap=30.)

For a conducting (or semimetallic materials whose electrons can be excited above the bandgap by the used pulse), you also need to define thermal parameters for the electronic system and their interaction with the phonons via electron-phonon-counpling:

my_conductor = SimMaterials(name='Blitzydummy', tdeb=300, cp_max=2.5e6, kappap=20., ce_gamma=100, kappae=100., gep=1e18)

For a magnetic material, whose spin dynamics you want to model with the M3TM, you need to define additional parameters within the model:

my_magnet = SimMaterials(name='Spinnydummy', tdeb=200, cp_max=2e6, kappap=10., ce_gamma=75, kappae=150., gep=0.8e18, spin=2.5, vat=1e-28, tc=600., muat=5., asf=0.06)

Create a sample structure

How to build a sample

With the materials you created before you can now build a sample.

my_sample = SimSample()

So far, our sampleholder is completely empty. Let's quickly grow a sample with a 5 nm Blitzydummy cap layer, a 15 nm Spinnydummy magnet and a 300 nm Gummydummy substrate.

You need to watch out for three things:

1. If you interface two materials, you always need to spcify boundary conditions for the interfacial phononic heat diffusion.
2. If you interface two materials with itinerant electrons, you also need to specify boundary conditions for the interfacial electronic heat diffusion.
3. You need two specify either a penetration depth for the laser pulse or a complex refractive index to simulate the penetration of the pump pulse into your sample structure(I will add both because I can!). For insulating materials introduce a penetration depth of exactly an integer 1!

Here is a list of the parameters to chose when adding layers to your sample:

• material (object). A material previously defined with the materials class
• dz (float). Layer thickness of the material in m. Important only for resolution of heat diffusion
• layers (int). Number of layers with depth material.dz to be added to the sample
• kappap_int [see 1.] (float/string). Phononic interface heat conductivity to the last block of the sample. Either in W/m/K or 'av', 'min', 'max' of the constants of the two interfaced materials
• kappae_int [see 2.] (float/string). Electronic interface heat conductivity to the last block of the sample. Either in W/m/K or 'av', 'min', 'max' of the constants of the two interfaced materials
• pen_dep [see 3.] (float). Penetration depth of the laser pulse in m if to be computed with Lambert-Beer absorption profile
• n_comp[see 3.] (complex float). Complex refractive index of the material. Use syntax 'n_r'+'n_i'j to initiate

my_sample.add_layers(material=my_conductor, dz= 1e-9, layers=5, pen_dep=7.5e-9, n_comp=2.8+8.5j)
my_sample.add_layers(material=my_magnet, dz=1e-9, layers=15, pen_dep=34e-9, n_comp=2.2+2.7j, kappae_int='max', kappap_int=1.)
my_sample.add_layers(material=my_insulator, dz=1e-9, layers=300, pen_dep=1, n_comp=1.2+0j, kappap_int='av')

Create a Pulse

How to create a Pulse

With the sample created we can now compute already how a pump pulse interacts with it.

To define the pulse you can/must introduce the following parameters:

• sample (object). Sample in use
• pulse_width (float). Sigma of gaussian pulse shape in s
• fluence (float). Fluence of the laser pulse in mJ/cm^2
• delay (float). Time-delay of the pulse peak after simulation start in s
• method (String). Method to calculate the pulse excitation map. Either 'LB' for Lambert-Beer or 'Abeles' for the matrix formulation calculating the profile via the Fresnel equations
• photon_energy_ev [only_Abeles] (float). Energy of the optical laser pulse in eV. Only necessary for method 'Abeles'
• theta [only_Abeles] (float). Angle of incidence of the pump pulse in respect to the sample plane normal in units of pi, so between 0 and 1/2. Only necessary for method 'Abeles'
• phi [only_Abeles] (float). Angle of polarized E-field of optical pulse in respect to incidence plane in units of pi, so between 0 and 1/2. Only necessary for method 'Abeles'

my_pulse_Abeles = SimPulse(sample=my_sample, method='Abeles', pulse_width=20e-15, fluence=5., delay=0.5e-12, photon_energy_ev=1.55, theta=1/4, phi=1/3)
my_pulse_Lambert_Beer = SimPulse(sample=my_sample, method='LB', pulse_width=20e-15, fluence=5., delay=0.5e-12)

You can then visualize the absorption profile with

my_pulse_Abeles.visualize(axis='z')

There are further options to visualize and save these figures, for documentation I refer to the documentation in the source code. Now, our pulse calculated with Abeles' matrix method looks like this:

Define Simulation Parameters

How to set up the simulation

All the physical parameters are defined now, we just need to run the simulation now. Therefor we need to define some computational parameters. These may influence if the simulation runs fast, slow, wonky or exact. In the end, you might need to play around a little until you have found a configuration that yields solid and fast results.

The parameters to define are:

• sample (object). The sample in use
• pulse (object). The pulse excitation in use
• end_time (float). Final time of simulation (including pulse delay) in s
• ini_temp (float). Initial temperature of electron and phonon baths in the whole sample in K
• solver (String). The solver used to evaluate the differential equation. See documentation of scipy.integrate.solve_ivp
• max_step (float). Maximum step size in s of the solver for the whole simulation
• atol [optional] (float). Absolute tolerance of solve_ivp solver. Default is 1e-6 as the default of the solver
• rtol [optional] (float). Relative tolerance of solve_ivp solver. Default is 1e-3 as the default of the solver

my_simulation = SimDynamics(sample=my_sample, pulse=my_pulse_Abeles, ini_temp=300., end_time=20e-12, solver='Radau', max_step=1e-13)

Let's just run it and see if we did okay here. To look at our data later we also need to save it on the harddrive. Automatically, in the package's regisrtry, a folder 'Results' will be created, where the simulation output files will be stored within a folder denoted by save_file. The actual data is stored in .npy format:

my_results = my_simulation.get_t_m_maps()
my_simulation.save_data(my_results, save_file='my_result_files')

Run the simulation

To run the simulation you have set up, open the terminal and navigate to ML_M3TM_Sim. Run the command

py -m code.Scripts.tutorial_input

The command prompt will show the setup of your sample, the computed pulse profile, initial conditions and briefly inform about the state of the simulation.

Sample constructed.
Constitiuents:  ['Blitzydummy' 'Spinnydummy' 'Gummydummy']
Thicknesses: [  5.  15. 300.]  nm
Number of layers: [5, 15, 300]

Absorption profile computed with Abeles' matrix method.
F =  5.0  mJ/cm^2
F_a_sim = 1.468017037261054 mJ/cm^2
F_r = 3.1144679343569615 mJ/cm^2
F_t = 0.4369338236503731 mJ/cm^2
F_a = F - F_r - F_t= 1.448598241992665 mJ/cm^2
Relative error due to finite layer size:  1.32 %

Starting simulation

Equilibration phase done.
Equilibrium magnetization in magnetic layers: [0.88880614 0.88880614 0.88880614 0.88880614 0.88880614 0.88880614
 0.88880614 0.88880614 0.88880614 0.88880614 0.88880614 0.88880614
 0.88880614 0.88880614 0.88880614]
at initial temperature: [300. 300. 300. 300. 300. 300. 300. 300. 300. 300. 300. 300. 300. 300.
 300.] K

Starting main simulation loop. Calculating until  19.5  ps after pulse maximum.
Solver:  Radau
maximum time step: 1e-13  s

Simulation done. Time expired: 46.86850047111511 s

Check the output

Within the plot module there are several methods to visualize your output data. Here is a quick guide to look at your data in some ways: