Model Background

Model Background

This page summarizes the physics and the numerical strategy behind the toolbox. 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 (TTM), 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 (surface point) 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 energy deposition with absorbance A, fluence F, temporal pulse profile g(t), and optical absorption coefficient alpha
• adiabatic boundary conditions 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 is S(t) = A * F * g(t) / Leff.

The two-stage multi-pulse strategy

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

The toolbox splits each pulse period into two stages:

1. Stage 1, pulse and relaxation. The full coupled TTM is integrated with MATLAB's stiff solver ode15s on a fine grid, from just before the pulse through electron-lattice equilibration (a few picoseconds).
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 pulse period costs a few hundred fine steps plus a fixed number of diffusion steps, so the approach scales to thousands of pulses. The baseline 50-pulse 0D run completes in about 0.4 s.

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 (Tl > Te).

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, Te undershoots Tl. The effect requires spatial transport: a 0D model cannot show it.

Depth_Profile_Solver reports per-pulse inversion metrics (magnitude, onset, duration), and 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 (Surface_Point_Solver, Radial_Profile_Solver, Scanning_Beam_Solver) uses a single total lattice conductivity because electron transport stays inside the heated layer. The 1D family (Depth_Profile_Solver, Single_Pulse_Visualizer) splits conduction into electron and phonon channels and adds the optical absorption coefficient.

Property	Unit	W	Cu	Au (0D only)	Al
gamma (electron heat capacity coeff.)	J m^-3 K^-2	137.3	98	67	136
Cl (lattice heat capacity)	J m^-3 K^-1	2.54e6	3.45e6	2.49e6	2.42e6
G (electron-phonon coupling)	W m^-3 K^-1	1.65e17	0.90e17	1.40e16	2.40e17
k total, 0D family	W m^-1 K^-1	174	401	317	237
ke0, 1D family	W m^-1 K^-1	150	390	n/a	220
kl, 1D family	W m^-1 K^-1	24	11	n/a	17
alpha (optical absorption)	m^-1	5.88e7	7.09e7	n/a	1.22e8

For tungsten, the solvers additionally use a temperature-dependent hybrid thermal conductivity during inter-pulse diffusion. Any metal can be modeled by setting material = 'custom' and supplying the manual property fields listed in Solver Reference.

Default absorbance is 0.55, appropriate for tungsten around 1030 nm at normal incidence; adjust it for your wavelength, surface condition, and material.

Assumptions to keep in mind

• One-dimensional heat flow in the depth solvers: valid while the spot radius is much larger than the thermal diffusion length, which holds for typical ~100 um spots on metals over MHz-scale gaps.
• Constant lattice heat capacity, no melting or phase-change latent heat. Predictions above the melt temperature indicate melting rather than resolve it.
• Simplified ke = ke0 * Te/Tl electron conductivity, a common approximation below the Fermi temperature.
• Beer-Lambert absorption with a fixed absorbance; no plasma or dynamic reflectivity effects at ablative fluences.
• Adiabatic boundaries: no radiative or convective surface losses (usually negligible on these timescales).

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 this toolbox implements, with parameter choices and experimental context for tungsten.

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