Guiding heat in laser ablation of metals on ultrafast timescales: an adaptive modeling approach on aluminum

Using an optimal control hydrodynamic modeling approach and irradiation adaptive time-design, we indicate excitation channels maximizing heat load in laser ablated aluminum at low energy costs. The primary relaxation paths leading to an emerging plasma are particularly affected. With impulsive pulses on ps pedestals, thermodynamic trajectories are preferentially guided in ionized domains where variations in ionization degree occur. This impinges on the gas-transformation mechanisms and triggers a positive bremsstrahlung absorption feedback. The highest temperatures are thus obtained in the expanding ionized matter after a final impulsive excitation, as the electronic energy relaxes recombinatively. The drive relies on transitions to weakly coupled front plasmas at the critical optical density, favoring energy confinement with low mechanical work. Alternatively, robust collisional heating occurs in denser regions above the critical point. This impacts the nature, the excitation degree and the energy content of the ablated matter. Adaptive modeling can therefore provide optimal strategies with information on physical variables not readily accessible and, as experimentally confirmed, databases for pulse shapes with interest in remote spectroscopy, laser-induced matter transfer, laser material processing and development of secondary sources.


Introduction
Energy coupling is a key element of many laser-matter interaction fields, ranging from secondary radiation sources or extreme matter states to laser material processing and remote sensing. Within this frame-work, increased attention is given to ultrafast laser sources, with the prospect of increasing the precision of interaction. The latter applications mentioned above mostly require controlled excitation and material removal events in ablation ranges (several J cm −2 ), as well as evolution paths that deliver heat confinement or highly excited matter. Specifically the generation of critical states that can spontaneously decompose into ionized phases is examined in applications where the emission of the laser ablation products has to be particularly intense, e.g. remote laser breakdown spectroscopy or radiation and particle sources [1,2]. This usually imposes the use of high-energy doses reflected in the photon cost. These states can nonetheless be achieved upon moderate laser exposure by correctly balancing the energy supply or the conversion between thermal and mechanical loads. The purpose is to obtain maximum effects from minimal energy costs related to the laser interaction. However, despite the fact that general ablation scenarios are well comprehended, a complete understanding of the recipe on how to obtain improved laser coupling and maximal excitation with its particular fine-tune control factors requires further investigation. A fundamental question is how to optimize the interaction by identifying ideal laser pulses for desired evolution scenarios and irradiation results [3]. Specifically related to ultrafast laser exposure, where matter properties vary from solid to plasma phases on electronic (τ e ) or ionic relaxation (τ hydro ) scales, impacting particularly the optical properties, the challenge is to improve the energetic balance beyond the standard material response. This can be regulated by acting on absorption channels and transitions dynamics with a modulated energy feedthrough. Well-defined thermomechanical paths can be induced and ablation aspects can be mastered to a considerable degree if irradiation is judiciously adjusted on these timescales.
Time manipulation of laser pulses has already been proven to have strong potential for achieving new synergies between light and matter. Temporal field gradients can determine enhanced plasma coupling and improved shock design [4][5][6]. Light confinement below the diffraction limit and nanostructuring can occur by energetically controlling carrier populations 3 with chirped pulses [7]. Coherent phonon populations may impact heat transport under the influence of multipulse sequences in phase with the vibration beat [8]. In laser ablated silicon Stoian et al [9] and Dachraoui and Husinsky [10] indicated a manyfold increase in ion emission with minimal energy costs under tailored ultrafast irradiation. The improvement was related to two factors: fast onset of absorptive liquid phases on picosecond scales or controlled generation of carrier plasmas on sub-picosecond scales. A similar ion enhancement was noted for metals [11] that involves hydrodynamic material expansion on tens of ps. This suggests a general effect of energy confinement assisting the phase transitions. Fast succession of density phases related to rarefaction was further evoked as a possible control source for nanodroplets in the ablation products [11,12]. This was related to Rayleigh-Taylor instabilities driven by the lifetime of the expelled liquid [12,13]. Equally, pulse forms were used to influence the excitation degree in laser-produced plasmas, either designed in simple double-pulse sequences or complex shapes determined by spectrally controlled feedback loops [12,14,15]. Significant spectral enhancement was thus obtained [12] relying on light exposure during the plasma formation stage and coupling in the incipient plasma phase. Alongside ablation products, pulse sequences could influence deformation rates down to almost complete evanescence of the rarefaction waves with minimal spallation [16,17]. In all these cases, laser-induced heat load and the thermomechanical balance play a fundamental role. The consequences of its eventual control are important for pressure-induced phase transitions, ablation characteristics and material removal efficiency, and we explore here global control options from a theoretical perspective.
In practice, experimental approaches offer limited access to primary evolution steps. In these conditions, a predictive theoretical method can provide extended information on how to achieve optimal interactions. Nonetheless, the complexity and lack of coherence of intermediate processes make it difficult to predict accurately essential improvement factors except by trials and educated guesses. Optimal control theories were then developed for controlling complicated parameter landscapes with light and have already been applied for regulating complex systems with coherent radiation [18][19][20]. Assuming the ansatz of controllability, i.e. the probability that a control target can be achieved with a laser controller [21], this allows, on the one hand, calculating optimal energy delivery [22,23] in the time domain for given objectives and, on the other hand, investigating active control knobs from temporal electric field distributions.
We extend the control approach to ablation-specific non-coherent phenomena occurring on picosecond scales and propose an adaptive hydrodynamic method by integrating numerical feedback loops into a hydrodynamic simulation code. With a focus on the temperature state parameter and its control, the procedure offers a predictive capability concerning thermodynamic evolutions with regulation laws not experimentally accessible. It will be shown that a numerical approach is powerful in predicting initial excitation conditions 4 required to increase and maximize the heat load in the ablation products, with major consequences for the ejected material states: excited species. The interest is twofold: an inquiry into regulation mechanisms and practical relevance for controlled ablation and excitation products. We are particularly interested here in ways of guiding laser energy towards obtaining higher energetic plasma phases in moderate irradiation regimes, beyond the levels obtainable by standard ultrafast laser ablation. A recipe for obtaining optimal pulse (OP) forms will be given and the 4 underlying physics will be discussed. This is motivated by the potential impact on the nature and on the energy content of the ablation species, with consequences for related applications.
The paper is organized as follows. The computational method is described in section 2, where details of the hydrodynamic model are given and the numerical feedback loop is described. The discussion section starts by evaluating the results of standard short pulse (SP) irradiation in terms of the limits of the achievable temperature levels. The results of the adaptive runs based on pulse time-design are then introduced, indicating the OP forms maximizing the temperature of the ablated matter in a significant range. This is accompanied by a discussion of the main mechanisms concurring with the optimal results and of the possible follow-ups in terms of applied and fundamental knowledge. Particular thermodynamic evolutions of ionized or dense plasma layers are followed, emphasizing heating scenarios in weakly coupled plasmas or supercritical fluids. The consequences of energetic states are pointed out. The conclusion section concludes the paper.

Hydrodynamic model
We modeled the nonequilibrium heating and expansion of aluminum (Al) under SP laser irradiation using a one-dimensional (1D) two-temperature non-equilibrium hydrodynamic code (Esther) [12,[24][25][26]. The use of a 1D model is justified for our intensity ranges and for spatial and temporal dimensions smaller than those associated with a Rayleigh-Taylor instability (µm, ns). As detailed hereafter, the Lagrangian approach solves the fluid equations for mass, momentum and energy conservation for electronic and ionic species upon energy absorption. Thus it allows a hydrodynamic and thermodynamic view of the evolution of the excited material. The phenomenological treatment is presented in detail in [25,26] and we only point out here the major points. The main steps of the hydrodynamic code are given below.
The interaction between the polarized pulsed laser field (λ = 800 nm) and the metal target is described by the Helmholtz wave equation in inhomogeneous media. The incident temporal laser electric fieldẼ(z = ∞, t) = |Ẽ inc (t)|e iω 0 t is injected as a boundary condition in the equation at normal incidence, being propagated towards negative z values: (1) Here˜ is the inhomogeneous complex dielectric permittivity of the metal, calculated as a function of electronic temperature T e , the ionic temperature T i , the density ρ and the average ionization Z . The electric field envelope can take various forms as will be seen in the next section.
Let us first discuss the optical response of the metal via its dielectric function. In considering absorption and the energy deposition rate, Al is viewed as a model case with parallel bands and a quasi-free-electron behavior in the solid phase. However, a deviation from a freeelectron-like behavior occurs around 1.5 eV: our irradiation conditions. The solid metal-light interaction includes therefore contributions from intraband (D, Drude-like free-electron heating) and interband (IB, Ashcroft-Sturm model [27]) electronic transitions, followed by free-electron absorption in the rarefied phases. The transient dielectric function˜ and the optical conductivitỹ σ are defined as where ω p = n e e 2 /m e 0 is the plasma frequency and ν eff is the effective (non-equilibrium) collision frequency. The latter is given by the sum of the electron-electron [ν ee (T e , n e )] and electron-phonon/ion [ν eph (T e , T i , n e )] [30] contributions in the solid phase and by [ν ei (T i , T e , n e )] in the plasma phase. We typically consider here the plasma phase as an ionized high-temperature low-density gas phase where the ionization is higher than 0.1. In calculating ν eff in the solid phase, the ν ee contribution originates from (e-e) umklapp collisional processes [25]. The electron-phonon contribution to scattering ν eph involves (e-ph) interactions in situations where the electronic temperature is larger than the ionic one. ν eph starts from a phenomenological collision frequency ν 0 eph which is fitted to match simultaneously tabulated data of n and k given by [28] and we add evolution laws depending on T e and T i [25]. The (e-ph) estimation follows the first-principle calculations from Lin et al [29] and the formalism of Kaganov et al [30]. In the plasma phase, ν ei = ν eq ei (T i , n e ) + ν neq ei (T e , n e ), indicating that ν ei scattering frequency is composed of an equilibrium and a non-equilibrium term. The equilibrium part ν eq ei (T i , n e ) is determined from the tabulated values of conductivity at equilibrium [31]. These tabulated values were shown to match well experimental situations dealing with detecting optical conductivities [32]. The second term corresponds to the nonequilibrium contribution to the collision frequency. For ν neq ei (T e , n e ), a Spitzer-like dependence has been assumed where the Coulomb logarithm has been corrected for high-density plasmas by the approach given in [34]. To ensure a convergence towards tabulated equilibrium values, we subtract the equilibrium contribution from the non-equilibrium term. In this manner, the contribution of electrons out of equilibrium is treated as a correction of the known values of the equilibrium conductivity.
Subsequently, equation (1) including the optical transient response (2) is jointly solved with the thermodynamic equations for the interacting electron and ion subsystems written in the Lagrangian form. The code explicitly considers the time-dependent thermal non-equilibrium (equations (3a) and (3b)). The specific energy conservation equations that take into account variations of the time-dependent energy exchange quantities in the z-direction are given below: with the terms detailed below: Here ε and K are the specific energy and thermal conductivity for electrons (e) and ions (i), respectively, and m is the mass of the Lagrangian cell. Additionally, V = 1/ρ is the specific volume of the calculation domain (cell) and p e−i the electrons or ions pressure.

6
The intervening terms shown explicitly in equations (4a)-(4e) are discussed below. Acting as the input source, Q L (z, t) is the specific optical power density source generated by the inverse bremsstrahlung absorption of the pump laser pulse, before being dissipated in the material. The energy transport and redistribution accounts for several influences [12]: (I) losses by diffusion Q e,i th ∼ ∂ z T e,i , (II) electron-phonon/electron-ion (e-i) collisional energy coupling with its specific exchange relaxation rate Q e−i , (III) the consumed mechanical work rate Q e,i mech = − p e,i ∂ t V and (IV) recombination kinetics Q rec following variations in the ionization states initiated by successive ionization/relaxation events. The absorption, relaxation and transfer energy terms [12] are material state dependent as explained below. The optical conductivitỹ σ relies on the permittivity˜ with [σ ] ∼ m [˜ ] determining the laser absorption.Ẽ represents the local field in interaction with a varying electronic population n e = n i Z . This path makes up an important energy exchange channel. As the electron density number varies, the extra energy is either supplied from the laser source (intrinsically considered via the optical conductivity) or released to the surrounding. This way it contributes to heating, with Q rec (z, t) indicated above accounting for the rate of specific energy recombination.
One important aspect should now be noted. Being driven essentially by temperature, these interaction terms affect the energy deposition and the transformation speed as several processes have antagonistic effects (i.e. temperature-driven collisional absorption or temperaturedependent effectiveness of heat diffusion). Just as an example, a femtosecond pulse exposure leads to high electronic temperatures that increase heat diffusion, in spite of the short duration of irradiation. This shows the importance of temperature control in the overall interaction scenario and, at the same time, suggests a possible optimization potential.
Two quantities have a determining role in establishing the thermal conversion and transport: the collisional e-ph/e-i coupling efficiency and a time-varying ionization state, both being strongly dependent on temperature and density. During evolution, the e-ph coupling initially determined from first-principle calculations [33] departs from the solid value [12] as γ (T e , n e ) = γ 0 (T e ) n e /n 0 e 1/3 where '0' denotes the solid phase (γ 0 (300 K) = 2.45 × 10 17 W m −3 K −1 ) and the symbols are usual quantities. For the dilute, hot and therefore weakly coupled plasmas, the e-i relaxation relies on the Spitzer model as discussed in [34]. This formalism accounts then for the relaxation-specific energy exchanged by e-ph and e-i collision Q e−i (z, t). At the same time, the average ionization Z varies during the electron heating phase consuming photons upon increase and releasing the energy by non-radiative recombination ∂ t Z | τ r during the quasi-adiabatic expansion. We recall that this impacts strongly on the absorption mechanisms and therefore on the energy storage. Assuming a strong offset from local equilibrium, the electronic energy drives the ionization values influencing thus the Q rec (z, t).
The Z values as a function of density (ρ) and electronic temperature (T e ) are taken from [31], being smeared by a characteristic decay time (τ r ) where a lower limit of 1 ps was imposed [35]. The energy that is delivered to the particle system as Z goes down corresponds to a variation in the electronic thermal energy and is reflected in the ionic temperature. The post-recombination states correspond to bound states in thermodynamic equilibrium with the surroundings. The model does not take into account the energy transfer and redistribution on specific atomic excitation levels, nor their particular bound-bound radiative relaxation, considered minimal in the irradiation conditions mentioned above. With the above dependences, we note that the macro-evolution of the ablation process is dominantly influenced by the specific energy content and finally by the ionic temperature (T i ), an essential factor that drives and, at the same time, is directly influenced by the succession of laser-induced effects. During relaxation, the material thermodynamic properties are described by an upgraded version of the Bushman-Lomonosov-Fortov multiphase equation of state spanning a large range of densities and temperatures between condensed states and hot plasmas. The values were reconstructed by Commissariatà l'Energie Atomique (CEA), France, based on [36]. Surface motion, induced shock wave and plasma expansion are considered by the hydrodynamic part of the calculation, satisfying the continuity relations given below: where u is the fluid velocity.

Numerical adaptive loop
To derive the best pulse form for correlated matter-light systems, a non-deterministic evolutionary search strategy [18] was developed. This relies on spectral phase modulation for tailoring the pulse shape [37]. An adaptive loop based on spectral phase modulation of an initial pulse E (ω) e iϕ(ω) with its corresponding time envelope E (t) was inserted in the hydrodynamic code. This is schematically shown in figure 1. A Fourier transfer function in the time domain is used to design the corresponding pulse forms for the code. Typically, the temporal envelopeẼ inc (t) of the laser electric field is obtained from the inverse Fourier transform of a Gaussian spectrum (full-width half-maximum (FWHM) λ = 2π c ω 2 ω = 6.2 nm, corresponding to a Fourier limited pulse duration of 150 fs), centered at the laser frequency ω 0 (recall that λ inc = 800 nm) and modulated by the spectral phase ϕ(ω) [37]: 8 Normalization to the laser fluence F is performed in order to ensure that, energetically, the same irradiation conditions are delivered by each tailored pulse, with c being the speed of light and 0 the dielectric vacuum permittivity.
The adaptive phase manipulation technique is briefly indicated below. The Fourier domain phase ϕ (ω) encompassing the spectral content of the pulse is considered as a string of 640 discrete values (pixels) with a spectral resolution of δλ = 0.05 nm. These values match typical values of standard commercial spatial light modulators (SLM). Being linked to the chosen spectral resolution, the shaping window extends to 42 ps, corresponding to pulses that can be designed in the laboratory with current optical modulators.
As mentioned, an evolutionary adaptive loop was generated between the programmable pulses and the hydrodynamic code, to improve the laser-generated outputs. As we are interested in the specific energy content, the desired optimal state is defined using a single physically relevant thermodynamic parameter, the ionic temperature T i , where experimental determination remains a challenge. The optimization algorithm is applied to examine the search space of the cost function f cost = max[T i (z, t)]. The optimality criterion, the maximum T i (irrespective of time and excitation zone) achieved in an irradiation sequence, serves therefore as a fitness parameter for the optimization procedure. Consequently, the phase is manipulated iteratively using genetic propagators, supporting the code with corresponding temporal forms that gradually impacts the temperature level. Via a pixel binning method, increasing time domains were explored, from less than 1 ps up to the maximal value. The convergence of the numerical procedure forecasts irradiation conditions for achieving the hottest states at constant energy requirements, allowing for prediction of the adequate control law. The result of the simulation does not only depend on the final state of the system but is also recursively determined by the history of the evolution paths during the regulated energy transport.

Standard short pulse (femtosecond) irradiation
Before approaching the optimization search, one has to establish the limits of the standard SP excitation. A first step in profiling the search space was made by exploring standard heating scenarios induced by femtosecond laser pulses (150 fs) using the hydrodynamic code. Typical temperatures achieved in different fluence ranges are thus calculated. Figure 2 shows the evolution of the maximal temperature T MAX of the hottest layers (HL) for a standard SP of 150 fs with increasing laser fluence F, mainly in ablation domains. Above 0.34 J cm −2 (the vaporization threshold at normal pressure) the material follows an ablative behavior with a sub-linear temperature increase. The value reached for a fluence of 3 J cm −2 , well in the ablation range, is approximately T SP = 22 × 10 3 K. The limited increase, localized at the expansion front, is particularly linked to the sudden SP energy deposition. Here the energy relaxation occurs primarily before the hydrodynamic expansion, at relatively constant quasisolid relaxation parameters, notably high specific heat that moderates the temperature gain. It has already been shown [12,38] that the strong initial compression supported by SP determines a swift mechanical expansion into the two-phase region, consuming the thermal energy in both expansion and the nucleation of the vapor phase.

Optimization results
The analysis of standard irradiation puts forward the question of whether, for the same energy input, the temperature can be further controlled. Under these circumstances the adaptive procedure was applied to find the extreme values that can be achieved by T i (z, t) when the pulse form and the energy delivery rate vary in time. As the temperature shows both time and space variations, we concentrated on the hottest zones and tried to influence the local maximal T i values in HL, identifying the highest and the least efficient pulses. An incident fluence of 3 J cm −2 in the ablation range was chosen. The results of the adaptive run show notably that, using the same energy input, it is possible to reach a far wider range of maximal temperatures, between 18 × 10 3 K (the lowest maximal value, denoted sometimes as minimal) and 55 × 10 3 K (the highest maximal value). The extended range of the temperature variation and the present discussion is specific to ablative regimes although observable effects in changing the highest temperature were obtained at lower fluences, below the threshold for achieving gasphase states. The pulses that deliver the low or high temperatures close to the extremes of the temperature ranges at this input fluence are depicted in figures 3(a) and (b). They indicate either a ps envelope with an impulsive ending for the highest T MAX or a short sequence for the lowest T MAX levels, respectively. We note therefore that an OP design (depicted in figure 3(a)) can determine a 2.5 times increase in temperature (see the OP level in figure 2). This relies on an energy delivery timescale up to the limit of the shaping window, in conditions where the incident fluence has not changed. Although the results are obtained at 3 J cm −2 input fluence, the observed pulse spread is also preserved in irradiation regimes beyond the chosen range to higher fluences, with apparently similar effectiveness (tested up to 10 J cm −2 ). The particular timescale suggests a departure from stress confinement conditions as the scale is limited by τ hydro λ D /c s , approximately 10 ps. Here λ D is the characteristic electronic diffusion length and c s is the speed of sound. This high temperature would require 10 times more energy with the SP at similar densities and shows 35% more efficiency than typical double-pulse sequences [12,15,16] or long Gaussian envelopes [12] within the same time frame. The lowest maximal temperature for the excited material is reached with an impulsive excitation, similar to SP, accompanied by some low-intensity wings ( figure 3(b)). This suggests a similar evolution to the one discussed in the first paragraphs, in conditions of shock confinement and mechanically driven expansion.
To capture the spatio-temporal dynamics, the evolution of the temperature in the expanding layers is given in figures 3(c) and (d) for the highest (MAX) and the least (MIN) efficient pulses. These graphs depict the levels of temperature achieved in the excited/ablated matter as the material expands in space and time.
We discuss below the case of the highest temperature. For the high temperature case the evolution of the highest temperature layer HL is indicated in figure 3(c). This shows that the maximum temperature is confined to the front of expansion. The highest value is obtained at the end of the sequence, a few ps after the maximum e-i non-equilibrium (marked as NE in figure 3(c)). As Al does not show strong optical variations as it makes the transition to the liquid phase [39], this points to an alteration of the properties related to the onset of the vapor phase. The gas-phase involvement is supported by the space evolution of the expanding matter on hundreds of nm, and it will be confirmed in the next paragraphs. We state, however, that the region of the highest temperature is not equivalent to the highest internal volume energy density (ρε i , marked as HE V ) located deeper towards the condensed phase. Apart from HL, initially placed at a depth of 2 nm (the initial position of the corresponding Lagrangian cell), a second expanding layer trajectory (10 nm) is given (SL) that illustrates a passage in the vicinity of the energy concentration region, reaching 42 × 10 3 K. The discussion concentrates first on HL as it locks the highest temperature. The particular relevance of the SL will become apparent later in the text (sections 3.3.2 and 3.3.4).

Optimal pulse characteristics: hydrodynamic scenario
The fundamental question is related to the specific OP time features that enable a maximal temperature, namely the low-intensity preliminary phase and the final impulse. The OP sequence shows a long pedestal of noisy low-intensity envelope for the first 41 ps. The apparent complicated energy spread can be approximated by a smooth constant supply of energy, as typically the sequence of low-intensity spikes are smoothed by a convolution with the ps material response. With an integral absorption twice as efficient as SP, this part of the pulse (80% of the total energy) is sufficient to raise the temperature up to 70% of the final OP value. A sudden well-pronounced peak is observed at the end of the sequence that achieves a final increase. Different optimization runs allow us to determine similar solutions where the obtained maximum temperature depends slightly on the number of iterations, within 2%. This suggests parameter space topologies with large maxima. The chosen solution has, nevertheless, the major characteristics of the ideal pulse.
We note that the optimal solution here was obtained with a constraint on the shaping window (limited to 42 ps). This is related to the practical feasibility of pulse shapes in current SLMs based on spectral filtering, with sub-nm spectral resolution. Improved solutions can be obtained on time scales extending to 150-200 ps, relying roughly on similar phenomena. In the next paragraphs, we indicate the particular states that enable this strong heating from a dynamic and energy balance standpoint.
A glance at the hydrodynamic progress of the excited material is given in figure 4. In parallel with the OP form ( figure 4(a)) superimposed on the instant absorbed power map (Q L ), it describes the sequence of density states ( figure 4(b)) and indicates for specific trajectories (HL, SL) the behavior of the speed of sound (figure 4(c)) and ionization degree (figure 4(d)), as indicators for the material state. The corroborated temperature and density information allows us to decompose the evolution into several stages.

3.3.1.
Gas-phase evolution. The first stage (0-10 ps) prepares the hydrodynamic advance to the gas phase (figures 4(a) and (b)). About 20% of the sequence energy is consumed to initiate vaporization (0.6 J cm −2 , approximately two times the threshold), equivalent to an absorbed fluence of 78 mJ cm −2 . The energy is deposited mainly in condensed phase over a depth √ 2Dτ = 28 nm, above the optical skin depth. This is sufficient to achieve an average specific energy in the range of 10 7 J kg −1 , energetically enough to initiate the plasma phase with a lowdensity evolution above the binodal (critical point (CP) energy density E C = 1.1 × 10 7 J kg −1 ). The density information ( figure 4(b)) at the end of this first stage of the laser exposure shows now the onset of several notable regions that continue to develop at later times: a low-density front and an increasingly dense material towards the depth. The vapor layer z sc , visible at about 10 ps and affecting the first nm of material, is optically transparent at 800 nm, with a subcritical electron density (n e < n c 10 21 cm −3 , z c being the critical density border at 800 nm). With its buffer function between the ablation products and the vacuum, it has the role of slowing down the evolution of the subsequent layers. As the expansion is repressed and the neighboring matter compressed, the underlying layers z int (among them HL) experience moderate heating for another 10 ps under the incoming energy and preserve a density high enough to allow nonnegligible e-i coupling. The HL optical density remains always superior to the critical value.
To understand what happens specifically in this time domain, we examine the dynamics of the phase transition by monitoring the material state-dependent speed of sound c s = √ (∂ p/∂ρ) S (figure 4(c)). This physical quantity reflects the ion pressure and serves as a sufficiently accurate probe for the eventual phase transitions. For HL c s stays roughly constant in the first picoseconds and then it experiences a drastic decrease at 9 ps as it penetrates into the two-phase region [40] towards vapor transformation. We note then a swift rise about 20 ps from a value of c s ∼ = 0 characteristic of the (practically non-ionized) liquid-vapor mixed phase in the vicinity of the CP to values typical of ionized gases, suggesting a rapid transition to an excited plasma phase characterized by low heat capacity and high optical conductivity. A similar dependence has the gas-phase expansion velocity mirroring the mechanical work. The HL enters the plasma region at the earliest time with an overall energy consumption in the liquid-vapor domain of 23% (the energy content of the pulse envelope between 9 and 20 ps corresponding to the liquid-vapor passage) under quasi-constant delivery that prevents cooling and allows an electronic temperature gain. During this, the liquid-vapor characteristic internal energy/particle augments due to e-i coupling, while the mechanical work stays low, comparable to thermal energy. As HL goes rapidly into a low-density high-temperature plasma phase with Z >1 (20 ps), this phase at n e n c further absorbs energy at double pace as light delivery continues between 20 and 42 ps. This stage involves 57% of the input sequence and accounts for 85% 13 of the total absorbed energy in the layer. With temperature, Z increases to quasi-solid values (figure 4(d)) in spite of the now accelerated expansion. The plasma is at present weakly coupled, with a descending scattering time. The bremsstrahlung absorption rises, acting in turn on energy deposition, increasing the accumulation. A positive feedback effect occurs reflected on the ultimate, maximal value of the ionic temperature around the end of the laser sequence. Here the effect of the final peak is deciding, releasing an impulsive heating of the plasma phase at the moment where absorbtion peaks (integral absorption at this moment is 68%). A short excitation is required to counteract the accelerated expansion of the hot plasma (still two times lower than in the case of SP excitation). Even if conductivity is high in this region, promoted by e-i nonequilibrium, the sudden ionization energy relaxes to the atomic system at the end of the sequence via recombination as Z decreases, being the primary factor for the temperature increase. This is equally supported by a moderate value of the heat capacity. This thermodynamic path facilitates then transient changes in the energy transport with the possibility of attaining states at T > T C .

Evolution of dense hot regions.
We have indicated before a second typical layer, SL, which advances not in the high-temperature zone but in the region of highest energy concentration HE V . As this maximizes the volume energy density accumulation which can be a reservoir for further evolution, its behavior is of interest and will be discussed below. The evolution of the SL and its adjacent regions follows a different path as compared to HL. Being localized deeper in supercritical electron density regions, where the laser light is screened (the electric field is a few times smaller (three times on average), as screening shows a progressive evolution from the solid phase (E HL /E SL 2) to the plasma phase (E HL /E SL 4)), the intermediate electronic heating stage is reduced. This happens despite a favorable hydrodynamic regime for increasing T related to a higher density. However, the energy density storage remains significant, particularly due to the high density. With restricted expansion due to z int , the SL does not show a typical gas-phase transition by crossing phase coexistence interfaces, as c s preserves a finite value, but evolves towards the plasma phase as a decomposing hot fluid with supercritical temperatures (see figure 3(c), T C = 6300 K for Al). The SL denotes therefore a typical supercritical behavior. Ionization Z remains limited, less sensitive to non-equilibrium, and the e-i collisions are the main factor for heating into warm dense states via supercritical paths.

Pulse analysis.
To investigate the importance of the presence of several peaks in the pulse form, we designed a specific method for analyzing OP, emphasizing the two types of energy-maximization behaviors discussed above: via thermodynamic gas-phase transitions and via supercritical paths. We defined a rectangular gate of approximately 2.5 ps which was scanned along the OP sequence replacing at each moment the local peaks with an equal but constant envelope energy domain. The results, given in figure 4(e), indicate that the small-amplitude peaks have little influence, mainly in keeping short-living non-equilibrium states. Nevertheless, the absence of the final peak impacting non-equilibrium T e is lowering the T MAX corresponding to the hottest HL temperature range, as this was related to plasma states at highest Z , even though the energy content of the laser sequence was not changed. This shows that the final contribution is related to a variation of Z and points out the importance of recombination in heating the HL, proceeding faster than the low-density collisional coupling. The SL is relatively insensitive to the final peak, confirming the alternative collisional heating mechanism. We have therefore indicated that the high-temperature regions in the proximity of HL and the high-energy-density SL regions have different driving elements. A summary of the relative importance of the various heating mechanisms (recombination or collisional heating) is given in figure 5, validating the scenario proposed above. This synthesizes the accumulation of specific energy via the different paths discussed so far, the ionization/recombination path ( t 0 Q rec (t ) dt ) and the electron-ion collisional path ( t 0 Q e−i (t ) dt ). In the superficial but optically dense HL layers ( figure 5(a)), the collisional heating mechanism is essential in the short-living solid phase, where the electronic energy relaxes via e-ph coupling in the first picoseconds after the interaction. The high electronic density favors a stronger collisional relaxation than recombination up to approximately 10 ps. The tendency is reversed as the material density drops during the liquid-gas transformation. The two contributions have similar weights in the liquid-vapor phase and at the beginning of the plasma phase. Here it is interesting to note that, between 10 and 25 ps, Q e−i and Q rec have similar magnitudes in equation (3a). This time window corresponds to relaxation rates of the order of about 10 ps as ∂ t Z /Z ∼ γ /c e where c e is the electronic specific heat capacity. An increase of the efficiency of the recombination path for heating the ions becomes visible after 25 ps for the OP. For the profound SL layers at higher densities ( figure 5(b)), despite the screened radiation, globally the collisional heating is most important as the accumulated effect shows at all times. This tendency is not affected by particular time moments where, for short instants, the instantaneous rate of recombination can become quite important, particularly at the last peak (not shown). It becomes evident that the dominant contribution for HL comes therefore from variation of the ionization degree and recombination, while the denser states of SL favor collisional transfer.
3.3.4. Thermodynamic phase space. We have discussed above variations of local parameters that we try now to position in the phase space in order to establish a global evolution. The global OP progress is summarized in the thermodynamic pressure-density-temperature p(ρ, T ) phase space ( figure 6(a)). The initial HL evolution deviates slightly from an isochoric transition due to its position close to the surface, which facilitates advancement without strong impediments. In addition, the pressure wave generated by the OP imposes a stress load distributed over the whole sequence (not shown), that does not force strong expansion up to the binodal crossing. Then, HL, although showing an evolution in the liquid-vapor zone, is able to rapidly leave the twophase zone re-intersecting the binodal at 20 ps and reaching a critical low-density plasma phase (ρ < ρ C ). The evolution crosses diagonally the space above the binodal for the last exposure part (half of the exposure time) into ionized absorbing states. Noteworthy is the SL behavior. The SL circumvents the CP after half of the irradiation sequence as high-pressure (GPa) supercritical fluid (justifying the acronym SL for a supercritical layer) but lacks the momentum for a strong evolution into the plasma phase. Following OP exposure 80% of the ablated matter follows supercritical paths as compared to SP (20%). A relevant indication of the OP supercritical paths (labeled as the initial positions of the Lagrangian cells) is given in figure 6(b), showing that, in addition to maximal heat loads, the OP equally favors supercritical paths as efficient ways towards the plasma phase.

Consequences of optimal coupling
The consequences of creating critical states are multiple, particularly in generating highly excited environments as requested e.g. in laser breakdown spectroscopy but also in secondary source development, accessing at these temperature values in EUV ranges. An example is offered in figure 7 showing the emitted integrated spectral energy density of the electrons within the first 100 ps. For this particular example the calculations were performed post-optimization by including a radiative transfer module. The integrated emissivity is determined by radiation transport. The equation of the radiative transfer in an expanding plasma is calculated with the assumption of an instantaneous distribution of emission and absorption sources (spatial zones that are radiating at different temperatures) determined by the local plasma conditions at each fluid element. A multigroup radiation transfer algorithm has been used with 100 energy groups of tabulated local thermodynamic equilibrium opacities extrapolated between the values of [28,41]. A spectral blueshift of 5 eV and a peak increase of more than one order of magnitude are to be seen for the OP. The optimal path found here could be general, as it relies on the achievement of high-conductivity plasma phases or supercritical expansion rather than on particular material properties. Essentially the pulse shows a material expansion and heating phase that is being driven into ionized domains where the rest of the energy is being deposited. This scenario also indicates a certain robustness against uncertainties in the code parameters. If the temperature level is sensitive to the chosen parameters which are known with some accuracy, the OP form drives essentially the same phenomena and preserves similar characteristics. Experimental confirmation of the benefits of longer, ps laser pulses determined by adaptive techniques based on either time-of-flight mass spectrometry or spectral emission have already been reported, confirming the effect of efficient plasma excitation [11][12][13] and the possibility of achieving high emissivity or spectral emission shifted towards higher energies. In these experiments, the monitored effect was a global (temperature-and densitydependent) plasma excitation behavior, spatially and temporally averaged for several hundreds of ns after the excitation, as probing temperature on the short time and scales reported here remains a challenge. However, they indicate as well the increased light coupling related to plasma interaction, and they validate thus the present results with the increased heating effect of longer envelopes on the emerging plasma. Secondly, as experimentally shown, favoring heat with respect to mechanical load and distributing it over a longer scale reduces tensile stress and cavitation. This impacts the sizes and density of nanodroplets in laser deposition techniques [12,13]. Thus, complementing the experimental approaches, the numerical optimal procedure supports a deeper understanding of the excitation control and access to primary relaxation steps. Finer excitation features are revealed, as no experimental noise or other limitations are to be considered. Depending on the chosen application, more developed fitness functions can be chosen to optimize the complex mix of several parameters. As a plus, the numerical adaptive loop allows us to address local parameters and reveals the fine details of excitation. It equally permits introspection into evolution paths non-accessible in experiments, giving a complementary view of the initial relaxation stages.

Conclusions
In conclusion, using a numerical adaptive approach, we have indicated that phase transitions can be guided by suitable energy delivery to maximum heat loads. The method focuses on primary relaxation events that establish further evolution in optimal conditions. Hot thermodynamic states are reached with minimal energy expenses by creating spatio-temporal correlations between flux delivery, absorption transients and relaxation, resulting from ps pedestals with impulsive endings. A plasma-based recipe was given which relies on early achievement of thermodynamically critical plasma phases at optical density. The role of recombination was outlined as a material-heating mechanism. Enhanced collisionally heated supercritical states, rich in volume energy, follow in denser regions, with a complementary collisional energy deposition mechanism. Adaptive modeling enables optimal irradiation strategies for applications in material processing and spectroscopy with new evolution paths and the prospect of reconfigurable irradiation tools.