Model Background

Model Background

This page summarizes the physics and the numerical strategy. The full treatment, including validation against experiments on tungsten, is in the companion article, doi:10.1007/s11665-026-14738-6.

The two-temperature model

On femtosecond timescales a laser pulse heats the electron gas of a metal long before the lattice responds. The two-temperature model, introduced by Anisimov and co-workers, tracks the electron and lattice baths separately and couples them through an electron-phonon coupling constant G.

The 0D form, solved by surface_point_solver and radial_profile_solver:

Ce(Te) * dTe/dt = -G*(Te - Tl) + S(t)      [electrons]
Cl     * dTl/dt = +G*(Te - Tl)             [lattice]

The 1D depth-resolved form, solved by depth_profile_solver and single_pulse_visualizer:

Ce(Te) * dTe/dt = d/dz[ ke(Te,Tl) * dTe/dz ] - G*(Te - Tl) + S(z,t)
Cl     * dTl/dt = d/dz[ kl * dTl/dz ]        + G*(Te - Tl)

with:

• Ce(Te) = gamma * Te, the Sommerfeld electronic heat capacity
• ke(Te,Tl) = ke0 * Te / Tl, a simplified Drude form of the electron thermal conductivity
• S(z,t) = A * F * g(t) * alpha * exp(-alpha*z), Beer-Lambert deposition with absorbance A, fluence F, temporal profile g(t), and optical absorption coefficient alpha
• adiabatic boundaries at the surface and the back face

In the 0D form the absorbed energy deposits into an effective heated thickness Leff, so the source term becomes S(t) = A * F * g(t) / Leff. When the post-pulse state is handed to the depth diffusion grid, the deposit is scaled so the grid receives exactly the layer's energy whatever its spacing. See the Validation page for how this differs from the MATLAB reference and how to reproduce the reference exactly.

Leff and the optical penetration depth are different quantities and are not interchangeable. Leff is a lumped thickness for the 0D energy balance and defaults to 100 nm. The depth-resolved solvers instead resolve absorption along depth through alpha_opt, whose reciprocal delta_opt is 17 nm for tungsten. Setting one when you meant the other is a common source of confusion, which is why the solvers accept only the key that applies to them.

The two-stage multi-pulse strategy

Resolving femtosecond pulse dynamics while simulating milliseconds of accumulation is the central difficulty of multi-pulse TTM. Brute force is prohibitive, because the stiff coupled equations need sub-picosecond resolution and a millisecond of laser time contains billions of such steps.

Each pulse period is split into two stages:

1. Stage 1, pulse and relaxation. The full coupled model is integrated on a fine grid from just before the pulse through electron-lattice equilibration, a few picoseconds. The depth-resolved solvers use a stiff BDF integrator, and the 0D family uses an adaptive RK4 scheme with a stability-limited step.
2. Stage 2, inter-pulse diffusion. Once the baths have equilibrated, the remaining gap to the next pulse is ordinary heat conduction. The end state of Stage 1 is mapped onto a coarser, deeper grid and advanced with an unconditionally stable Crank-Nicolson scheme.

The surface temperature at the end of Stage 2 becomes the initial condition of the next pulse. Each period costs a few hundred fine steps plus a fixed number of diffusion steps, so the approach scales to thousands of pulses. The 50-pulse 0D baseline completes in well under a tenth of a second.

The surface temperature inversion

A signature prediction of the depth-resolved model is a transient inversion at the surface, where the lattice becomes hotter than the electrons.

After the pulse, surface electrons are extremely hot, and the large electron thermal conductivity drains their energy into the bulk quickly. The lattice meanwhile heats through electron-phonon coupling at rate G. When electron heat leaves the surface faster than coupling replaces it, the electron temperature undershoots the lattice temperature. The effect requires spatial transport, so a 0D model cannot show it.

depth_profile_solver reports per-pulse inversion metrics covering magnitude, onset, and duration. inversion_quantifier aggregates them across a multi-pulse run, including how the inversion evolves as base temperature accumulates.

Material presets

The solvers ship with literature-based presets. The 0D family, meaning the surface point, radial profile, and scanning beam solvers, uses a single total conductivity, because electron transport stays inside the heated layer. The depth-resolved family splits conduction into electron and phonon channels and adds the optical absorption coefficient.

Property	Unit	W	Cu	Al	Au
gamma, electron heat capacity coefficient	J m^-3 K^-2	137.3	98	136	67
Cl, lattice heat capacity	J m^-3 K^-1	2.54e6	3.45e6	2.42e6	2.49e6
G, electron-phonon coupling	W m^-3 K^-1	1.65e17	0.90e17	2.40e17	1.40e16
total conductivity, 0D family	W m^-1 K^-1	174	401	237	317
ke0, depth-resolved family	W m^-1 K^-1	150	390	220	not available
kl, depth-resolved family	W m^-1 K^-1	24	11	17	not available
alpha_opt, optical absorption	m^-1	5.88e7	7.09e7	1.22e8	not available
melting point	degC	3422	1085	660	1064

The total conductivity of the 0D family is ke0 + kl, which is why the two shapes carry the same physics.

Gold has no published electron and phonon conductivity split or absorption coefficient in this table, so the depth-resolved solvers reject it rather than guess. It works in the 0D family.

For tungsten the solvers additionally use a measured temperature-dependent conductivity curve during inter-pulse diffusion, running from 174.9 W/mK at 293 K down to 95.0 at 2500 K. Supplying your own conductivity replaces that curve with a constant value, since the curve would otherwise silently ignore what you asked for. Set kTable to measured or constant to control that choice explicitly.

Any metal can be modelled by setting material to custom and supplying the manual property keys listed in Solver Reference. To change a single property of an existing preset, pass that property directly instead, for example {"material": "W", "delta_opt": 23e-9}.

Default absorbance is 0.55, appropriate for tungsten near 1030 nm at normal incidence. Adjust it for your wavelength, angle, surface condition, and temperature.

Assumptions to keep in mind

• One-dimensional heat flow in the depth solvers. This holds while the spot radius is much larger than the thermal diffusion length, which is the usual case for 100 micron spots on metals over MHz-scale gaps. The depth solver reports its lateral diffusion length and warns when that assumption weakens.
• Constant lattice heat capacity, with no melting and no latent heat. A prediction above the melting point tells you the parameter set melts the surface, it does not resolve the melt.
• The simplified ke = ke0 * Te/Tl electron conductivity, a common approximation below the Fermi temperature.
• Beer-Lambert absorption with a fixed absorbance, with no plasma or dynamic reflectivity at ablative fluences.
• Adiabatic boundaries, with no radiative or convective surface losses. These are negligible on the timescales modelled here.

References

• S. I. Anisimov, B. L. Kapeliovich, T. L. Perel'man, "Electron emission from metal surfaces exposed to ultrashort laser pulses," Sov. Phys. JETP 39, 375 (1974). The original two-temperature model.
• D. Fieser, U. N. Dewanjee, A. Hu, "A Computationally Efficient Two-Stage Two-Temperature Model for Multi-pulse Femtosecond Laser Heat Accumulation in Tungsten," Journal of Materials Engineering and Performance (2026). doi:10.1007/s11665-026-14738-6. The model implemented here, with parameter choices and experimental context for tungsten.

See the article's bibliography for the wider TTM literature.