Landau-Lifshitz-Bloch simulations#
Here we calculate the vectorial magnetization dynamic in a magnetic heterostructure employing the mean-field quantum Landau-Lifshitz-Bloch (LLB) apporach. Please read the following review to get an overview of the LLB equation in ultrafast magnetism
U. Atxitia, D. Hinzke, and U. Nowak,
Fundamentals and Applications of the Landau-Lifshitz-Bloch Equation,
J. Phys. D. Appl. Phys. 50, (2017).
Here we need to solve the following differential equation:
The three terms describe
precession at Larmor frequency,
transversal damping (conserving the macrospin length), and
longitudinal damping (changing macrospin length due to incoherent atomistic spin excitations within the layer the macrospin is defined on).
\(\alpha_{\parallel}\) and \(\alpha_{\perp}\) are the longitudinal damping and transverse damping parameters, respectively. \(\gamma_e = -1.761\times10^{11}\,\mathrm{rad\,s^{-1}\,T^{-1}}\) is the gyromagnetic ratio of an electron.
The effective magnetic field is the sum of all relevant magnetic interactions:
where
\(\mathbf{H}_\mathrm{ext}\) is the external magnetic field
\(\mathbf{H}_\mathrm{A}\) is the uniaxial anisotropy field
\(\mathbf{H}_\mathrm{ex}\) is the exchange field
\(\mathbf{H}_\mathrm{th}\) is the thermal field
The definitions of these subterms are described in their respective API documentation.
The material parameters of the current example arbitrarily chosen.
Setup#
Do all necessary imports and settings.
import udkm1Dsim as ud
u = ud.u # import the pint unit registry from udkm1Dsim
import scipy.constants as constants
import numpy as np
import time
import matplotlib.pyplot as plt
%matplotlib inline
u.setup_matplotlib() # use matplotlib with pint units
Structure#
Refer to the structure-example for more details.
We are providing a magnetization of the atoms which are used as initial magnetization in the LLB simulations, in case no specific initial condition is provided by the user.
Co = ud.Atom('Co', mag_amplitude=1, mag_gamma=0*u.deg, mag_phi=0*u.deg)
Ni = ud.Atom('Ni', mag_amplitude=1, mag_gamma=90*u.deg, mag_phi=0*u.deg)
Fe = ud.Atom('Fe', mag_amplitude=1, mag_gamma=0*u.deg, mag_phi=90*u.deg)
Si = ud.Atom('Si')
Solving the LLB requires several additional parameters of the according Layer
objects:
eff_spin
- effective spincurie_temp
- Curie temperaturelamda
- parameter used in the longitudinal & transversal damping terms \(\alpha_\parallel\) & \(\alpha_\perp\) (misspelled because of pythonlambda
-functions)mag_moment
- atomic magnetic momentaniso_exponent
- exponent of the uniaxial anisotropyanisotropy
- vector of the anisotropy \([K_x, K_y, K_z]\)exch_stiffness
- vector of the exchange stiffness \([A_{i\rightarrow (i-1)}, A_{i\rightarrow i}, A_{i\rightarrow (i+1)}]\)mag_saturation
- zero temperature saturation magnetization
The correct physical units and more details are documented in the Layer API.
prop_Ni = {}
# two-temperture model
prop_Ni['heat_capacity'] = ['0.1*T',
532*u.J/u.kg/u.K,
]
prop_Ni['therm_cond'] = [20*u.W/(u.m*u.K),
80*u.W/(u.m*u.K),]
g = 4.0e18 # electron-phonon coupling
prop_Ni['sub_system_coupling'] = \
['-{:f}*(T_0-T_1)'.format(g),
'{:f}*(T_0-T_1)'.format(g)
]
prop_Ni['lin_therm_exp'] = [0, 11.8e-6]
prop_Ni['opt_ref_index'] = 2.9174+3.3545j
# LLB parameters
prop_Ni['eff_spin'] = 0.5
prop_Ni['curie_temp'] = 630*u.K
prop_Ni['lamda'] = 0.005
prop_Ni['mag_moment'] = 0.393*u.bohr_magneton
prop_Ni['aniso_exponent'] = 3
prop_Ni['anisotropy'] = [0.45e6, 0.45e6, 0.45e6]*u.J/u.m**3
prop_Ni['exch_stiffness'] = [5e-14, 5e-14, 5e-14]*u.J/u.m
prop_Ni['mag_saturation'] = 500e3*u.J/u.T/u.m**3
# build the layer
layer_Ni = ud.AmorphousLayer('Ni', 'Ni amorphous', thickness=1*u.nm,
density=7000*u.kg/u.m**3, atom=Ni, **prop_Ni)
Number of subsystems changed from 1 to 2.
# similar to Ni layer
prop_Co = {}
prop_Co['heat_capacity'] = ['0.1*T',
332*u.J/u.kg/u.K,
]
prop_Co['therm_cond'] = [20*u.W/(u.m*u.K),
80*u.W/(u.m*u.K),]
g = 5.0e18
prop_Co['sub_system_coupling'] = \
['-{:f}*(T_0-T_1)'.format(g),
'{:f}*(T_0-T_1)'.format(g)
]
prop_Co['lin_therm_exp'] = [0, 11.8e-6]
prop_Co['opt_ref_index'] = 2.9174+3.3545j
prop_Co['eff_spin'] = 3
prop_Co['curie_temp'] = 1480*u.K
prop_Co['lamda'] = 0.005
prop_Co['mag_moment'] = 0.393*u.bohr_magneton
prop_Co['aniso_exponent'] = 3
prop_Co['anisotropy'] = [0.45e6, 0.45e6, 0.45e6]*u.J/u.m**3
prop_Co['exch_stiffness'] = [5e-14, 5e-14, 5e-14]*u.J/u.m
prop_Co['mag_saturation'] = 1400e3*u.J/u.T/u.m**3
layer_Co = ud.AmorphousLayer('Co', 'Co amorphous', thickness=1*u.nm,
density=7000*u.kg/u.m**3, atom=Co, **prop_Co)
Number of subsystems changed from 1 to 2.
# similar to Ni layer
prop_Fe = {}
prop_Fe['heat_capacity'] = ['0.1*T',
732*u.J/u.kg/u.K,
]
prop_Fe['therm_cond'] = [20*u.W/(u.m*u.K),
80*u.W/(u.m*u.K),]
g = 6.0e18
prop_Fe['sub_system_coupling'] = \
['-{:f}*(T_0-T_1)'.format(g),
'{:f}*(T_0-T_1)'.format(g)
]
prop_Fe['lin_therm_exp'] = [0, 11.8e-6]
prop_Fe['opt_ref_index'] = 2.9174+3.3545j
prop_Fe['eff_spin'] = 2
prop_Fe['curie_temp'] = 1024*u.K
prop_Fe['lamda'] = 0.005
prop_Fe['mag_moment'] = 2.2*u.bohr_magneton
prop_Fe['aniso_exponent'] = 3
prop_Fe['anisotropy'] = [0.45e6, 0.45e6, 0.45e6]*u.J/u.m**3
prop_Fe['exch_stiffness'] = [5e-14, 5e-14, 5e-14]*u.J/u.m
prop_Fe['mag_saturation'] = 200e3*u.J/u.T/u.m**3
layer_Fe = ud.AmorphousLayer('Fe', 'Fe amorphous', thickness=1*u.nm,
density=7000*u.kg/u.m**3, atom=Fe, **prop_Fe)
Number of subsystems changed from 1 to 2.
# this is the non-magnetic substrate
prop_Si = {}
prop_Si['heat_capacity'] = [100*u.J/u.kg/u.K, 603*u.J/u.kg/u.K]
prop_Si['therm_cond'] = [0, 100*u.W/(u.m*u.K)]
prop_Si['sub_system_coupling'] = [0, 0]
prop_Si['lin_therm_exp'] = [0, 2.6e-6]
prop_Si['sound_vel'] = 8.433*u.nm/u.ps
prop_Si['opt_ref_index'] = 3.6941+0.0065435j
layer_Si = ud.AmorphousLayer('Si', "Si amorphous", thickness=1*u.nm, density=2336*u.kg/u.m**3,
atom=Si, **prop_Si)
Number of subsystems changed from 1 to 2.
S = ud.Structure('NiCoFeNi')
S.add_sub_structure(layer_Ni, 10)
S.add_sub_structure(layer_Co, 15)
S.add_sub_structure(layer_Fe, 10)
S.add_sub_structure(layer_Ni, 25)
S.add_sub_structure(layer_Si, 200)
S.visualize()

Initialize Heat and the Excitation#
h = ud.Heat(S, True)
h.save_data = False
h.disp_messages = True
h.excitation = {'fluence': [25]*u.mJ/u.cm**2,
'delay_pump': [0]*u.ps,
'pulse_width': [0.15]*u.ps,
'multilayer_absorption': True,
'wavelength': 800*u.nm,
'theta': 45*u.deg}
# temporal and spatial grid
delays = np.r_[-10:10:0.05, 10:200:0.5]*u.ps
_, _, distances = S.get_distances_of_layers()
Calculate Heat Diffusion for 2-Temperature Model#
# enable heat diffusion
h.heat_diffusion = True
# set the boundary conditions
h.boundary_conditions = {'top_type': 'isolator', 'bottom_type': 'isolator'}
# The resulting temperature profile is calculated in one line:
temp_map, delta_temp = h.get_temp_map(delays, 300)
Surface incidence fluence scaled by factor 0.7071 due to incidence angle theta=45.00 deg
Calculating _heat_diffusion_ for excitation 1:1 ...
Absorption profile is calculated by multilayer formalism.
Total reflectivity of 46.5 % and transmission of 4.0 %.
Elapsed time for _heat_diffusion_ with 1 excitation(s): 20.389669 s
Calculating _heat_diffusion_ without excitation...
Elapsed time for _heat_diffusion_: 35.628406 s
Elapsed time for _temp_map_: 56.187758 s
plt.figure(figsize=[6, 8])
plt.subplot(2, 1, 1)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, temp_map[:, :, 0],
shading='auto')
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('Temperature Map Electrons')
plt.subplot(2, 1, 2)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, temp_map[:, :, 1],
shading='auto')
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('Temperature Map Phonons')
plt.tight_layout()
plt.show()

Landau-Lifshitz-Bloch Simulations#
The LLB
class requires a Structure
object and a boolean force_recalc
in order to overwrite previous simulation results.
These results are saved in the cache_dir
when save_data
is enabled.
Printing simulation messages can be en-/disabled using disp_messages
and progress bars can be enabled using the boolean switch progress_bar
.
llb = ud.LLB(S, True)
llb.save_data = False
llb.disp_messages = True
print(llb)
Landau-Lifshitz-Bloch Magnetization Dynamics simulation properties:
Magnetization simulation properties:
This is the current structure for the simulations:
Structure properties:
Name : NiCoFeNi
Thickness : 260.00 nanometer
Roughness : 0.00 nanometer
----
10 times Ni amorphous: 10.00 nanometer
15 times Co amorphous: 15.00 nanometer
10 times Fe amorphous: 10.00 nanometer
25 times Ni amorphous: 25.00 nanometer
200 times Si amorphous: 200.00 nanometer
----
no substrate
Display properties:
================ =======
parameter value
================ =======
force recalc True
cache directory ./
display messages True
save data False
progress bar True
================ =======
Brillouin Function#
Internally, the LLB calculates a mean-field magnetization map for the according electron temperatures
\(T_e\) at for every layer and for every time step. This is done by solving the Brillouin function of each
layer and then mapping the result onto the according spatio-temporal grid, as given by the temp_map
.
mean_field_mag_map = llb.get_mean_field_mag_map(temp_map[:, :, 0])
Calculating _mean_field_magnetization_map_ ...
Elapsed time for _mean_field_magnetization_map_: 1.736737 s
plt.figure(figsize=[6, 8])
plt.subplot(2, 1, 1)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, temp_map[:, :, 0],
shading='auto')
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('Temperature Map Electrons')
plt.subplot(2, 1, 2)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, mean_field_mag_map, shading='auto')
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('Mean-Field magnetization')
plt.tight_layout()
plt.show()

In order to run the actual LLB simulation, an optional initial magnetization init_mag
as well as an external magnetic field H_ext
can be provided.
While H_ext
is provided in cartesian coordinates \([H_x, H_y, H_z]\) in Tesla, the initial magnetization is given in polar coordinates \([A, \phi, \gamma]\) in units of [none, rad, rad], following the definitions shown in the user guide.
Running the simulation is done by calling the get_magnetization_map
method which requires a delay
vector as well as the temp_map
.
The resulting magnetization_map
is returned in polar coordinates.
init_mag = np.array([1.0, (0.*u.deg).to('rad').magnitude, (0*u.deg).to('rad').magnitude])
magnetization_map = llb.get_magnetization_map(delays, temp_map=temp_map, init_mag=init_mag,
H_ext=np.array([0, 0.05, 1]))
Calculating _magnetization_map_ ...
Calculating _mean_field_magnetization_map_ ...
Elapsed time for _mean_field_magnetization_map_: 1.756780 s
Elapsed time for _LLB_: 2.209331 s
Elapsed time for _magnetization_map_: 2.210315 s
plt.figure(figsize=[6, 12])
plt.subplot(3, 1, 1)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, magnetization_map[:, :, 0],
shading='auto', cmap='RdBu_r', vmin=-1, vmax=1)
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('$A$')
plt.subplot(3, 1, 2)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, magnetization_map[:, :, 1],
shading='auto', cmap='RdBu_r', vmin=-3.14, vmax=3.14)
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('$\phi$')
plt.subplot(3, 1, 3)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, magnetization_map[:, :, 2],
shading='auto', cmap='RdBu_r', vmin=-3.14, vmax=3.14)
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('$\gamma$')
plt.tight_layout()
plt.show()

plt.figure(figsize=[6,8])
plt.subplot(2,1,1)
plt.plot(delays, np.mean(magnetization_map[:, 0:10, 0], axis=1), label=r'$A$')
plt.legend()
plt.xlabel('Delay (ps)')
plt.ylabel('Magnetization')
plt.subplot(2,1,2)
plt.plot(delays, (np.mean(magnetization_map[:, 0:10, 1], axis=1)*u.rad).to('deg'), label=r'$\phi$')
plt.plot(delays, (np.mean(magnetization_map[:, 0:10, 2], axis=1)*u.rad).to('deg'), label=r'$\gamma$')
plt.legend()
plt.xlabel('Delay (ps)')
plt.ylabel('Magnetization')
plt.show()

The helper
modules provides functions to convert the result also into cartesian coordinates:
magnetization_map_xyz = ud.helpers.convert_polar_to_cartesian(magnetization_map)
plt.figure(figsize=[6, 12])
plt.subplot(3, 1, 1)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, magnetization_map_xyz[:, :, 0],
shading='auto', cmap='RdBu', vmin=-1, vmax=1)
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('$M_x$')
plt.subplot(3, 1, 2)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, magnetization_map_xyz[:, :, 1],
shading='auto', cmap='RdBu', vmin=-1, vmax=1)
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('$M_y$')
plt.subplot(3, 1, 3)
plt.pcolormesh(distances.to('nm').magnitude, delays.to('ps').magnitude, magnetization_map_xyz[:, :, 2],
shading='auto', cmap='RdBu', vmin=-1, vmax=1)
plt.colorbar()
plt.xlabel('Distance [nm]')
plt.ylabel('Delay [ps]')
plt.title('$M_z$')
plt.tight_layout()
plt.show()

plt.figure(figsize=[6,8])
plt.subplot(3,1,1)
plt.plot(delays, np.mean(magnetization_map_xyz[:, 0:10, 0], axis=1), label=r'$M_x$ Ni')
plt.plot(delays, np.mean(magnetization_map_xyz[:, 25:35, 0], axis=1), label=r'$M_x$ Fe')
plt.legend()
plt.xlabel('Delay (ps)')
plt.ylabel('Magnetization')
plt.subplot(3,1,2)
plt.plot(delays, np.mean(magnetization_map_xyz[:, 0:10, 1], axis=1), label=r'$M_y$ Ni')
plt.plot(delays, np.mean(magnetization_map_xyz[:, 25:35, 1], axis=1), label=r'$M_y$ Fe')
plt.legend()
plt.subplot(3,1,3)
plt.plot(delays, np.mean(magnetization_map_xyz[:, 0:10, 2], axis=1), label=r'$M_z$ Ni')
plt.plot(delays, np.mean(magnetization_map_xyz[:, 25:35, 2], axis=1), label=r'$M_z$ Fe')
plt.legend()
plt.xlabel('Delay (ps)')
plt.ylabel('Magnetization')
plt.show()

In both subplots it is obvious that the magnetization already starts to change before the actual laser excitation raises the electronic temperatures.
This is because of the misalignment of the initial magnetization and the external magnetic field.
To avoid this, the steady-state solution of the layer magnetization must be found first and used as init_mag
for the actual transient simulations.