Using experimental values of heat capacity with complicated temperature dependence

[mamattern]

I would like to use for the electronic heat capacity the sum of the magnetic and electronic contribution, which can be extracted from the experimentally determined total heat capacity. For modelling the transient temperatures for small (below the Curie temperature) and large fluences (electron temperature exceeds the Curie temperature) it is necessary to describe the heat capacity on both sides of the critical temperature. Instead of finding a function for the complete temperature range that reasonably describes the experimental heat capacity, I would like to directly use the experimental values for the heat capacity in the toolbox.

My Idea is using np.interp(T, temperature_list, hat_capacity_list) which returns the heat capacity for the temperature T. Following Issue #105 I try to define the heat caapcity and the integral by:

prop_Ni['heat_capacity'] = [1, 1]
Ni_layer = Ni = ud.UnitCell('Ni, 'Ni layer', caxis_Ni, prop_Ni)

Ni_layer._heat_capacity =  ['lambda T: 1', 'lambda T: 300']

my_temp = [300, 400, 500, 600, 700]
my_int_heat_capacity = [300, 400, 500, 600, 700]

Ni_layer.int_heat_capacity = [ 'interp(T, my_temp, my_int_heat_capacity)', '300*T']

Unfortunately this does not work and I receive an error message:

  File "udkm1dsim_ni20nm_2tm.py", line 173, in <module>
    temp_map, delta_temp_map = h.get_temp_map(delays, Init_temp)

  File "udkm1Dsim\simulations\heat.py", line 739, in get_temp_map
    self.calc_temp_map(delays, init_temp)

  File "udkm1Dsim\simulations\heat.py", line 870, in calc_temp_map
    temp, _ = self.get_temperature_after_delta_excitation(fluence,

  File "udkm1Dsim\simulations\heat.py", line 700, in get_temperature_after_delta_excitation
    final_temp[i, 0] = brentq(fun, init_temp[i, 0], 1e5)

  File "scipy\optimize\_zeros_py.py", line 783, in brentq
    r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)

  File "udkm1Dsim\simulations\heat.py", line 697, in fun
    return (masses[idx]*(int_heat_capacities[idx][0](final_temp)

  File "<lambdifygenerated-668>", line 2, in _lambdifygenerated
    [T_0, T_1] = _Dummy_116

TypeError: cannot unpack non-iterable float object

[dschick]

Hi Max,

when you set the private attributes starting with the _ you should directly use the lambda-functions and not their string representation.
Please also set the Ni_layer._int_heat_capacity directly to bypass the automatic lambda-fication.
This snippet should work

Ni_layer._heat_capacity = [lambda T: 1, lambda T: 300]
Ni_layer._int_heat_capacity = [lambda T: np.interp(T, my_temp, my_int_heat_capacity), lambda T: 300*T]

[mamattern]

Thank you for the answer and the advice for the correct syntax. Unfortunately, I get an error message again. But this time a different one:

  File "udkm1dsim_ni20nm_2tm.py", line 175, in <module>
    temp_map, delta_temp_map = h.get_temp_map(delays, Init_temp)

  File "udkm1Dsim\simulations\heat.py", line 739, in get_temp_map
    self.calc_temp_map(delays, init_temp)

  File "udkm1Dsim\simulations\heat.py", line 870, in calc_temp_map
    temp, _ = self.get_temperature_after_delta_excitation(fluence,

  File "udkm1Dsim\simulations\heat.py", line 700, in get_temperature_after_delta_excitation
    final_temp[i, 0] = brentq(fun, init_temp[i, 0], 1e5)

  File "scipy\optimize\_zeros_py.py", line 783, in brentq
    r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)

ValueError: f(a) and f(b) must have different signs

[dschick]

okay, this seems to be a more fundamental numerical issue.
I assume that you use 0 pulse duration in your simulations?
In that case the get_temperature_after_delta_excitation() method is called to calculate the temperature after delta excitation.
This numerical minimizer seem to have a problem, if the heat capacity is not monotonous.

Could you simply try using finite pulse durations? This should avoid calling this method.

[mamattern]

For finite pulse lengths, the code in your previous message works perfectly. Surprisingly, when I checked the results quantitatively, I found that giving the integral is not needed. In fact, the integral has no influence on the result. If I give an integral significantly higher than expected from the given heat capacity, the transient temperature does not change.

Is this expected from your point of view?

[dschick]

Indeed this behavior is more than expected.

If you call the get_temp_map without a pulse duration of the excitation (independent of heat diffusion enabled/disabled), the temperature jump at the arrival of the pump laser is always calculated by get_temperature_after_delta_excitation().
In order to calculate the temperature change of, e.g., an UnitCell due to the amount of absorbed energy, one needs to integrate its heat_capacity up to the final temperature:

This is actually the only reason why the int_heat_capacity is implemented.
If you are using finite pulse durations in combination with heat diffusion, the pump laser absorption is included via the source term in the NTM differential equations:

I am still wondering about your last error message:

ValueError: f(a) and f(b) must have different signs

To my understanding, this happens when the int_heat_capacity is non-monotonous.
If I am not mistaken your heat_capacity is non-monotonous because of the peak at the magnetic phase transition.
But the int_heat_capacity should be monotonous again!?
Could you maybe paste the according plots of both here?

I was also thinking if it is easier to replace the analytical integration with a numerical integration.
Obviously the analytical version is quite error prone and complex.

[mamattern]

Thank you for clarifying this difference, it was not on my radar.

Regarding the problem with the failing integration of the heat capacity in the case of 0 pulse duration:
This problem only occurs, if the integral does not increase with temperature with a value below the deposited energy. I think the problem is, that there is no solution of the equation determining the final temperature, if the deposited energy exceeds the maximum value of the given integral of the heat capacity. This happened in my first example, because I did not consider that np.interp(400, [100, 200, 300], [100, 200, 300]) returns 300 as integral for T=400, which was the case for my simple testing example.

In summary,

my_heat_capacity = np.genfromtxt('eq_model.dat')/58.69e-3
my_int_heat_capacity = np.genfromtxt('eq_model_integral.dat')/58.69e-3

Ni_layer._heat_capacity = [lambda T: np.interp(T, my_temp, my_heat_capacity), lambda T: 442.1]
Ni_layer._int_heat_capacity = [lambda T: np.interp(T, my_temp, my_int_heat_capacity), lambda T: 442.1*T]

works perfectly for both finite pulse durations and 0 pulse duration. Also for the heat capacity and its integral below. Thank you for the support.

[dschick]

cool that this is working now. I will close this issue but will think about removing the analytical integral at all, see #109