Method to set the initial magnetization to equilibrium value prior to excitation

[AvonReppert]

In the newly implemented LLB class that solves the Landau Lifshitz Bloch equation it would be very useful to have a method that determines the equilibrium magnetization direction for a given external field.

Using this method one could initialize the magnetization state to the equilibrium value prior to laser-excitation. One possibility would be to run the LLB simulation without laser-excitation with an increased damping until the change in the magnetization between two timesteps falls below a certain threshold or a time-limit is reached.

In the process of implementing the Landau-Lifshitz-Gilbert LLG equation that accounts also for the magneto-elastic coupling we run into the same issue, which we approach by letting the simulation run for many nanoseconds prior to the laser excitation at t = 0 so that the magnetization reaches has reached its equilibrium value as shown in the following image:

This process is however very costly as it creates huge magnetization map arrays with irrelevant data that only reflects the precession towards the equilibrium state starting from a user-defined definition. Since the same equilibration problem affects the LLB and the m3TM I wonder if you could provide a method that initializes the magnetization map of this structure to its equilibrium value for a given initial Temperature, initial strain state and initial external field direction.

[AvonReppert]

We need to allow the direct import of a given magnetization profile so we need to slightly modify the get_magnetization method of the llb class in the magnetization module. Currently it only allows for a vector with 3 entries but Nx3 entries where N is the number of layers should also be accepted.

            if 'init_mag' in kwargs:
                if not isinstance(kwargs['init_mag'], np.ndarray):
                    raise TypeError(
                        'init_mag must be a numpy ndarray with all elements representing the number of layers.')

                distances, _, _ = self.S.get_distances_of_layers(False)
                N = len(distances)
                init_mag_shape = kwargs['init_mag'].shape

                if init_mag_shape not in [(3,), (N, 3)]:
                    raise ValueError(
                        f'init_mag must be a vector with shape (3,) or (N, 3) with N being the number of layers.')

[friweber]

I solved the problem, at least for the case where only one magnetic layer exists. I did this by defining

def calc_initial_magnetization_map(self, delays_pre, layer_size, **kwargs):
        
            strain_pre = np.zeros((len(delays_pre),layer_size))
            temp_pre = np.zeros((len(delays_pre),layer_siz

            temp_pre += 300 #initial temperature 
        
        
            magnetization_map = self.calc_magnetization_map(delays_pre,temp_pre, strain_pre, **kwargs)
            
            
            t1 = time()
            
            self.disp_message('Elapsed time for initial_magnetization_map:'
                              ' {:f} s'.format(time()-t1))
            return magnetization_map[-1]

delays_pre is an array for defining the time before time zero and layer size is the number of unit cells in the single magnetic layer. Then strain_pre and temp_pre are initialized. By using calc_magnetization_map, the magnetization dynamic is calculated. The last entries in the calculated magnetization_map are returned, as this is the equilibrium magnetization. Here is an example how I used it:

init_mag = np.array([1, 0, 0])
  
init_mag_map = llg.calc_initial_magnetization_map(delays_pre, length, 
                                                  H_ext=np.array([0.4*np.sin(angle), 0., 0.4*np.cos(angle)]), init_mag=init_mag)

init_mag = np.array([np.mean(init_mag_map[:, 0]), np.mean(init_mag_map[:, 1]), np.mean(init_mag_map[:, 2])])

magnetization_map = llg.get_magnetization_map(delays, temp_map=temp_map, strain_map=strain_map, H_ext=np.array([
                                              0.4*np.sin(angle), 0., 0.4*np.cos(angle)]), init_mag=init_mag)

e))

As one can see, the calculated initial_magnetization is the equilibrium magnetization since there is nothing happening before time zero.

[AvonReppert]

Agree something along this line is what I have in mind.

The more I look into it, the harder it is to provide a generalized function that allows to reach the equilibrium state prior to excitation. Maybe it is best that we write a function that works in the scenarios we are using with the LLG equation in this branch and add it into this class. With the LLB model I still receive error messages when I try to insert an artificially generated homogeneous temp-map.

I was thinking to make a method out of this code snipped, which work for the LLG example but somehow fail for the LLB example:

init_temp = 300
delays_equilibrium = np.r_[-100:1000:1]*u.ps
temp_map_equilibrium = np.full((len(delays_equilibrium), len(distances), 2), init_temp)
H_ext = np.array([1, 0.05, 1])
is_magnetic = S.get_layer_property_vector('_curie_temp')>0

threshold = 2e-5
for counter in range(3):

    if counter == 0:
        init_mag = np.array([1.0, (0.*u.deg).to('rad').magnitude, (0*u.deg).to('rad').magnitude])
        magnetization_map_equilibrium = llb.get_magnetization_map(delays_equilibrium, temp_map=temp_map_equilibrium, init_mag=init_mag,
                                              H_ext=np.array([0, 0.05, 1]))
    else:
        init_mag = magnetization_map_equilibrium[-1, :, :]
        magnetization_map_equilibrium = llb.get_magnetization_map(delays_equilibrium, temp_map=temp_map_equilibrium, init_mag=init_mag,   H_ext=np.array([0, 0.05, 1]))

    deviations = np.sum(np.abs(magnetization_map_equilibrium[-1, :, :]-magnetization_map_equilibrium[-2, :, :]))
    print("residual deviation after " + str(counter) + " iterations: " + str(deviations))
    if np.sum(is_magnetic) > 0:
        mean_A = np.mean(magnetization_map_equilibrium[-1, is_magnetic, 0])
        mean_phi = np.mean(magnetization_map_equilibrium[-1, is_magnetic, 1])*u.rad
        mean_gamma = np.mean(magnetization_map_equilibrium[-1, is_magnetic, 2])*u.rad

        print("Amplitude = " + str(mean_A) + r" phi = " + str(mean_phi.to('deg')) + r" gamma = " + str(mean_gamma.to('deg')))
    else:
        print("The sample contains no magnetic layers")

    if deviations < threshold:
        break

init_mag = magnetization_map_equilibrium[-1, :, :]

magnetization_map = llb.get_magnetization_map(delays, temp_map=temp_map, init_mag=init_mag,
                                              H_ext=np.array([0, 0.05, 1]))

however, if I try it with the LLB I get errors that I am still trying to understand.

[friweber]

If I understand you correctly, your Idea is to implement something that checks whether the magnetization is stable, and then returns the magnetization. So it's not just calculating 10000ps but rather stops after i.e. 2000 since the magnetization is stable?

[AvonReppert]

Yes, I would like to have a method that runs the equilibration N times in steps of 1ps and at the end of each run checks if the magnetization is stable (i.e. all changes are below a user defined threshold). The user sets the threshold, and the number N of iterations until the magnetization is stable. Just to avoid that the precession doesn't cease after a sensible amount of time because the damping is too low. The method should print the coordinates of the average A, phi and gamma after each iteration just to see for the user if the magnetization is stable.

[dschick]

Hi @AvonReppert & @friweber

many thanks for your inputs!

Currently it only allows for a vector with 3 entries but Nx3 entries where N is the number of layers should also be accepted.

This is of course a first step and I have not thought about that in the beginning - but it should be very easy to implement.

however, if I try it with the LLB I get errors that I am still trying to understand.

Could you please post some of them in order to get some insights?

your Idea is to implement something that checks whether the magnetization is stable, and then returns the magnetization.

I guess we can do that already within the ode solver of python:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

There one can define an event after which the solver should terminate. This event could be reaching a certain stability condition. Hence one can simply run the simulation for a rather very long delay range and let the solver do the job.

At the same time one should check the atol and rtol parameters.

@Nilodirf also suggested temporally increasing the damping parameters of the sample stack to faster converge to the equilibrium state.

[AvonReppert]

Could you please post some of them in order to get some insights?

Well essentially I wonder why the following minimal "example" case does not run in the LLB class. I wanted to provide a homogenous temperature map initialized at 300K directly, but it does not work but rather runs into mathematical errors.

import matplotlib.pyplot as plt
import numpy as np
import udkm1Dsim as ud
u = ud.u  # import the pint unit registry from udkm1Dsim
u.setup_matplotlib()  # use matplotlib with pint units


# %% Layer properties

# Loading atoms
Ni = ud.Atom('Ni', mag_amplitude=1, mag_gamma=90*u.deg, mag_phi=0*u.deg)
Si = ud.Atom('Si')

prop_Ni = {}
# parameters for a 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.5
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'] = [0.1e-15, 1e-15, 0.01e-15]*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)


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)


# %% Building the structure

S = ud.Structure('Ni')
S.add_sub_structure(layer_Ni, 20)
S.add_sub_structure(layer_Si, 200)
S.visualize()

# %% Heat simulation without excitation
init_temp = 300  # K
fluence = 0     # mJ/cm^2
# temporal and spatial grid
delays = np.r_[-100:200:1]*u.ps
_, _, distances = S.get_distances_of_layers()


h = ud.Heat(S, force_recalc=False)
h.save_data = True
h.disp_messages = True
h.excitation = {'fluence': [fluence]*u.mJ/u.cm**2,
                'delay_pump':  [0]*u.ps,
                'pulse_width':  [0]*u.ps,
                'multilayer_absorption': True,
                'wavelength': 800*u.nm,
                'theta': 45*u.deg}

h.heat_diffusion = True
h.boundary_conditions = {'top_type': 'isolator', 'bottom_type': 'isolator'}
temp_map, delta_temp = h.get_temp_map(delays, init_temp)


# %% LLB simulation without  excitation

H_ext = np.array([1, 0, 1])
delays_equilibrium = np.r_[-100:200:1]*u.ps
temp_map_equilibrium = np.full((len(delays_equilibrium), len(distances), 2), init_temp)
is_magnetic = S.get_layer_property_vector('_curie_temp') > 0

llb = ud.LLB(S, True)

llb.save_data = True
llb.disp_messages = True

mean_field_mag_map_1 = llb.get_mean_field_mag_map(temp_map_equilibrium[:, :, 0])
init_mag_1 = np.array([1, (0.*u.deg).to('rad').magnitude, (0*u.deg).to('rad').magnitude])


magnetization_map_1 = llb.get_magnetization_map(
    delays_equilibrium, temp_map=temp_map_equilibrium, init_mag=init_mag_1, H_ext=H_ext)
magnetization_map_xyz_1 = ud.LLB.convert_polar_to_cartesian(magnetization_map_1)

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


# %% Reduce the damping by a factor 1000
prop_Ni['lamda'] = 0.0005
layer_Ni = ud.AmorphousLayer('Ni', 'Ni amorphous', thickness=1*u.nm,  density=7000*u.kg/u.m**3, atom=Ni, **prop_Ni)

S = ud.Structure('Ni')
S.add_sub_structure(layer_Ni, 20)
S.add_sub_structure(layer_Si, 200)

# %% Heat simulation with excitation
init_temp = 300  # K
fluence = 15     # mJ/cm^2
_, _, distances = S.get_distances_of_layers()


h = ud.Heat(S, force_recalc=False)
h.save_data = True
h.disp_messages = True
h.excitation = {'fluence': [fluence]*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}

# 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_2, delta_temp_2 = h.get_temp_map(delays, init_temp)


# %%
llb = ud.LLB(S, True)

llb.save_data = True
llb.disp_messages = True


mean_field_mag_map_2 = llb.get_mean_field_mag_map(temp_map_2[:, :, 0])
init_mag_2 = magnetization_map_1[-1, :]


magnetization_map_2 = llb.get_magnetization_map(delays, temp_map=temp_map_2, init_mag=init_mag_2, H_ext=H_ext)
magnetization_map_xyz_2 = ud.LLB.convert_polar_to_cartesian(magnetization_map_2)


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

It yields error messages like this:

udkm1Dsim\udkm1Dsim\simulations\magnetization.py:1049: RuntimeWarning: divide by zero encountered in divide
  dBdx = 1 / (two_eff_spins**2 * np.sinh(x / (two_eff_spins))**2) \
udkm1Dsim\udkm1Dsim\simulations\magnetization.py:1050: RuntimeWarning: divide by zero encountered in divide
  - (two_eff_spins + 1)**2 / \
udkm1Dsim\udkm1Dsim\simulations\magnetization.py:1049: RuntimeWarning: invalid value encountered in subtract
  dBdx = 1 / (two_eff_spins**2 * np.sinh(x / (two_eff_spins))**2) \
code\udkm1Dsim\udkm1Dsim\simulations\magnetization.py:1219: RuntimeWarning: invalid value encountered in cast
  chi_long[under_tc] = np.divide(
code\udkm1Dsim\udkm1Dsim\simulations\magnetization.py:970: RuntimeWarning: divide by zero encountered in divide
  H_th[under_tc] = 1/(2 * chi_long[under_tc]) * (
code\udkm1Dsim\udkm1Dsim\simulations\magnetization.py:971: RuntimeWarning: divide by zero encountered in divide
  1 - mag_map_squared[under_tc]/mf_magnetizations[under_tc]**2
code\udkm1Dsim\udkm1Dsim\simulations\magnetization.py:970: RuntimeWarning: invalid value encountered in cast
  H_th[under_tc] = 1/(2 * chi_long[under_tc]) * (
code\udkm1Dsim\udkm1Dsim\simulations\magnetization.py:1091: RuntimeWarning: divide by zero encountered in divide
  np.divide(lambdas[under_tc], mf_magnetizations[under_tc]), (
code\udkm1Dsim\udkm1Dsim\simulations\magnetization.py:1090: RuntimeWarning: invalid value encountered in cast
  alpha_trans[under_tc] = np.multiply(

The funny thing is, that a code where I generate the temp_map from an excitation with 0 fluence works although the resulting temp_map is identical. So the following works:

import matplotlib.pyplot as plt
import numpy as np
import udkm1Dsim as ud
u = ud.u  # import the pint unit registry from udkm1Dsim
u.setup_matplotlib()  # use matplotlib with pint units


# %% Layer properties

# Loading atoms
Ni = ud.Atom('Ni', mag_amplitude=1, mag_gamma=90*u.deg, mag_phi=0*u.deg)
Si = ud.Atom('Si')

prop_Ni = {}
# parameters for a 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.5
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'] = [0.1e-15, 1e-15, 0.01e-15]*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)


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)


# %% Building the structure

S = ud.Structure('Ni')
S.add_sub_structure(layer_Ni, 20)
S.add_sub_structure(layer_Si, 200)
S.visualize()

# %% Heat simulation without excitation
init_temp = 300  # K
fluence = 0     # mJ/cm^2
# temporal and spatial grid
delays = np.r_[-100:200:1]*u.ps
_, _, distances = S.get_distances_of_layers()


h = ud.Heat(S, force_recalc=False)
h.save_data = True
h.disp_messages = True
h.excitation = {'fluence': [fluence]*u.mJ/u.cm**2,
              'delay_pump':  [0]*u.ps,
              'pulse_width':  [0]*u.ps,
              'multilayer_absorption': True,
              'wavelength': 800*u.nm,
              'theta': 45*u.deg}

h.heat_diffusion = True
h.boundary_conditions = {'top_type': 'isolator', 'bottom_type': 'isolator'}
temp_map, delta_temp = h.get_temp_map(delays, init_temp)


# %% LLB simulation without  excitation

H_ext = np.array([1, 0, 1])
is_magnetic = S.get_layer_property_vector('_curie_temp') > 0

llb = ud.LLB(S, True)

llb.save_data = True
llb.disp_messages = True

mean_field_mag_map_1 = llb.get_mean_field_mag_map(temp_map[:, :, 0])
init_mag_1 = np.array([1, (0.*u.deg).to('rad').magnitude, (0*u.deg).to('rad').magnitude])


magnetization_map_1 = llb.get_magnetization_map(delays, temp_map=temp_map, init_mag=init_mag_1, H_ext=H_ext)
magnetization_map_xyz_1 = ud.LLB.convert_polar_to_cartesian(magnetization_map_1)

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


# %% Reduce the damping by a factor 1000
prop_Ni['lamda'] = 0.0005
layer_Ni = ud.AmorphousLayer('Ni', 'Ni amorphous', thickness=1*u.nm,  density=7000*u.kg/u.m**3, atom=Ni, **prop_Ni)
S = ud.Structure('Ni')
S.add_sub_structure(layer_Ni, 20)
S.add_sub_structure(layer_Si, 200)

# %% Heat simulation with excitation
fluence = 15     # mJ/cm^2
_, _, distances = S.get_distances_of_layers()


h = ud.Heat(S, force_recalc=False)
h.save_data = True
h.disp_messages = True
h.excitation = {'fluence': [fluence]*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}

# 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_2, delta_temp_2 = h.get_temp_map(delays, init_temp)


# %%
llb = ud.LLB(S, True)

llb.save_data = True
llb.disp_messages = True


mean_field_mag_map_2 = llb.get_mean_field_mag_map(temp_map_2[:, :, 0])
init_mag_2 = magnetization_map_1[-1, :]


magnetization_map_2 = llb.get_magnetization_map(delays, temp_map=temp_map_2, init_mag=init_mag_2, H_ext=H_ext)
magnetization_map_xyz_2 = ud.LLB.convert_polar_to_cartesian(magnetization_map_2)


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

and this I do not understand.

[friweber]

It has been some time since I wrote anything here, but I have made progress. I have built a function that finds the inertial magnetisation relatively quickly. In the following figure, you can see that before time 0, the magnetisation is constant, and after the excitation, the dynamics starts.

The function looks like this; for a period of 0-500 ps, the dynamics are calculated, if the magnetisation no longer changes, then the inertial magnetisation is reached, if the inertial magnetisation is not reached, then the calculation continues with the average value of the last calculated dynamics. This function is part of the magnetization class

def calc_equilibrium_magnetization_map(self, H_ext = np.array([0, 0, 0]), threshold = 1e-6, iterations = 100,  **kwargs):
        
        t1 = time()
        
        distances, _, _ = self.S.get_distances_of_layers(False)
        is_magnetic = self.S.get_layer_property_vector('_curie_temp')>0
        
        delays_equilibrium = np.linspace(0, 500, 1000)*u.ps
        delays_equilibrium = delays_equilibrium.to('s').magnitude
        
        strain_equilibrium = np.zeros((len(delays_equilibrium), len(distances)))
        
        M = len(delays_equilibrium)
        
        init_temp = 300.
        temp_equilibrium = np.full((len(delays_equilibrium), len(distances), 2), init_temp)
        
        
        
        for counter in range(iterations):

            if counter == 0:
                init_mag = H_ext/np.linalg.norm(H_ext)
                init_mag = self.convert_cartesian_to_polar(init_mag)
                magnetization_map_equilibrium = self.calc_magnetization_map(delays_equilibrium,temp_equilibrium, strain_equilibrium, H_ext, init_mag, **kwargs)
            else:
                init_mag_A = np.mean(magnetization_map_equilibrium[:, :, 0], axis = 0)
                init_mag_phi = np.mean(magnetization_map_equilibrium[:, :, 1], axis = 0)
                init_mag_gamma = np.mean(magnetization_map_equilibrium[:, :, 2], axis = 0)
                
                init_mag = np.zeros_like(magnetization_map_equilibrium[-1, :, :])
                init_mag[:,0] = init_mag_A 
                init_mag[:,1] = init_mag_phi
                init_mag[:,2] = init_mag_gamma 
                
                magnetization_map_equilibrium = self.calc_magnetization_map(delays_equilibrium,temp_equilibrium, strain_equilibrium, H_ext, init_mag, **kwargs)
        
            deviations_A = np.mean(magnetization_map_equilibrium[-1, is_magnetic, 0]) - np.mean(magnetization_map_equilibrium[-2, is_magnetic, 0])
            deviations_phi = np.mean(magnetization_map_equilibrium[-1, is_magnetic, 1]) - np.mean(magnetization_map_equilibrium[-2, is_magnetic, 1])
            deviations_gamma = np.mean(magnetization_map_equilibrium[-1, is_magnetic, 2]) - np.mean(magnetization_map_equilibrium[-2, is_magnetic, 2])
            deviations = np.sqrt(deviations_A**2+deviations_phi**2+deviations_gamma**2)
            print("residual deviation after " + str(counter) + " iterations: " + str(deviations))
            if np.sum(is_magnetic) > 0:
                mean_A = np.mean(magnetization_map_equilibrium[-1, is_magnetic, 0])
                mean_phi = np.mean(magnetization_map_equilibrium[-1, is_magnetic, 1])*u.rad
                mean_gamma = np.mean(magnetization_map_equilibrium[-1, is_magnetic, 2])*u.rad
        
                print("Amplitude = " + str(mean_A) + r" phi = " + str(mean_phi.to('deg')) + r" gamma = " + str(mean_gamma.to('deg')))
            else:
                print("The sample contains no magnetic layers")
        
            if deviations < threshold:
                break

        init_mag = magnetization_map_equilibrium[-1, :, :]

        
        self.disp_message('Elapsed time for initial_magnetization_map:'
                          ' {:f} s'.format(time()-t1))
        return init_mag

A new problem arises, however, which possibly has to do with the boundary conditions. I simulated 20nm nickel on 200nm silicon. It turned out that strange things happen at the edge between silicon and nickel.

The code for the simulation looks like this:

import matplotlib.pyplot as plt
import numpy as np
import udkm1Dsim as ud
u = ud.u  # import the pint unit registry from udkm1Dsim
u.setup_matplotlib()  # use matplotlib with pint units


# %% Layer properties

# Loading atoms
Ni = ud.Atom('Ni', mag_amplitude=1, mag_gamma=90*u.deg, mag_phi=0*u.deg)
Si = ud.Atom('Si')

prop_Ni = {}
# parameters for a 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.0005
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'] = [0.1e-15, 1e-15, 0.01e-15]*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)


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)


# %% Building the structure

S = ud.Structure('Ni')
S.add_sub_structure(layer_Ni, 20)
S.add_sub_structure(layer_Si, 200)
S.visualize()

# %% Heat simulation with excitation
init_temp = 300  # K
fluence = 15     # mJ/cm^2
_, _, distances = S.get_distances_of_layers()


h = ud.Heat(S, force_recalc=False)
h.save_data = True
h.disp_messages = True
delays = np.r_[-100:200:0.1]*u.ps
h.excitation = {'fluence': [fluence]*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}

# enable heat diffusion
h.heat_diffusion = True
# set the boundary conditions
h.boundary_conditions = {'top_type': 'isolator'}
# The resulting temperature profile is calculated in one line:

temp_map, delta_temp = h.get_temp_map(delays, init_temp)

pnum = ud.PhononNum(S, True)
strain_map = pnum.get_strain_map(delays, temp_map, delta_temp)

is_magnetic = S.get_layer_property_vector('_curie_temp')>0
# %%
llb = ud.LLB(S, True)

llb.save_data = True
llb.disp_messages = True

H_ext = np.array([1,0,1])


init_mag_map = llb.calc_equilibrium_magnetization_map(H_ext=H_ext, threshold = 1e-5)


magnetization_map = llb.get_magnetization_map(delays, temp_map=temp_map, strain_map = strain_map, init_mag=init_mag_map, H_ext=H_ext)
magnetization_map_xyz = ud.LLB.convert_polar_to_cartesian(magnetization_map)


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

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.xlim(0,20)
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.xlim(0,20)
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.xlim(0,20)
plt.ylabel('Delay [ps]')
plt.title('$M_z$')

plt.tight_layout()
plt.show()

[dschick]

Hi @friweber ,

Thanks for the continuation! Maybe you should move the second part of your example to a dedicated issue and if so, it would be good to know, whether the problem depends on the thickness of the magnetic layer.

Although your code regarding the equilibrium magnetization seems to work properly, I am wondering if the time steps/duration could be adapted somehow?
Have you also tried temporally increasing the damping of the magnetic layers?
This might help to find the equilibrium much faster.

Thanks

Daniel