# Physics-informed Machine Learning for Better Understanding Laser-Matter Interaction

Fayad Ali Banna, Jean-Philippe Colombier, Rémi Emonet, Marc Sebban

► To cite this version:

Fayad Ali Banna, Jean-Philippe Colombier, Rémi Emonet, Marc Sebban. Physics-informed Machine Learning for Better Understanding Laser-Matter Interaction. The 36th IEEE International Conference on Tools with Artificial Intelligence (ICTAI 2024), IEEE, Oct 2024, Herndon, VA, United States. pp.7. hal-04717837

HAL Id: hal-04717837

https://hal.science/hal-04717837v1

Submitted on 18 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Physics-informed Machine Learning for Better Understanding Laser-Matter Interaction

Fayad Ali Banna*†, Jean-Philippe Colombier*, Rémi Emonet‡*†, and Marc Sebban*‡

*Université Jean Monnet Saint-Etienne, CNRS, Institut d'Optique Graduate School, Laboratoire Hubert Curien UMR 5516, F-42023 St-Etienne, France

†Institut Universitaire de France, ‡Inria

**Abstract**—Physics-informed machine learning typically assumes that the underlying physical laws are known and abundant training data is available. These assumptions do not hold in the context of self-organization of matter, a phenomenon that leads to the emergence of patterns when a surface is irradiated with an ultrafast laser beam. Indeed, due to the constraints of the electronic data acquisition devices, the creation of large datasets is made impossible. Moreover, modeling this dynamic process is challenging as it involves coupling between electromagnetism, thermodynamics and fluid mechanics under far-from-equilibrium conditions that are not yet fully understood. This paper aims at taking a step forward towards a better understanding of this complex phenomenon. We specifically focus on the laser energy absorption of the surface, which is governed by the distinctive characteristics of Maxwell's equations in an inhomogeneous lossy medium. This involves modelling physics at the nano scale and incurs high simulation costs that make any exploration impractical. To address this major issue, we investigate different physics-informed learning models. In this low data regime, our study reveals that learning a simple U-Net-based surrogate model surpasses (i) more sophisticated neural architectures and (ii) the FDTD-based solver in speed by several orders of magnitude. Interestingly, our study highlights a link between the formation of patterns and the magnitude of absorbed energy.

**Index Terms**—Machine learning, Knowledge discovery, Physics-informed ML, Maxwell's equations.

## I. INTRODUCTION

In the real world, physical phenomena are often governed by complex systems described by partial differential equations (PDEs). Despite important scientific advances in numerical simulation, solving efficiently PDEs often remains prohibitively expensive. To address this issue, Physics-informed machine learning (PiML) [3], [4] has recently emerged as a compelling strategy for (i) crafting efficient surrogate solvers, (ii) augmenting physical laws with data-driven models or more generally (iii) constraining physics aware models [8]. Building upon this available physical knowledge, many algorithmic contributions have allowed to (i) build suitable physics-regularized loss functions, (ii) guide the initialization of deep neural networks with consistent parameters or (iii) design new

(data+knowledge)-driven hybrid neural architectures. This new line of research led to a novel generation of deep-learning frameworks aiming to predict either the solution of a PDE or the dynamics itself (e.g., PINN [12], FNO [5], PINO [6], U-NO [11], PDE-Net [9], to cite a few).

Current PiML methods often assume well-known underlying physical laws or the availability of abundant training data. The scenario, where (i) the PDEs only partially explain the observations and (ii) there is a scarcity of data, has been minimally explored up to now in the PiML literature. An example reflecting this unfavorable setting is self-organization of matter, the target application considered in this paper. Self-organization is the spontaneous process where some form arises from an initially disordered state. Nature provides many such examples (e.g., zebra stripes, ripples of sand dunes, sea shells, etc.). In surface engineering, self-organization of matter occurs when an ultra-short laser pulse interacts with a surface (called light coupling), inducing the emergence of patterns by material modification (see Fig. 1). Indeed, when the pulse interacts with the surface, the material heats up, becomes liquid and moves, leading to the appearance of regularity that can have significant practical value, e.g., in optics [16], microbiology [10], or to induce hydrophobicity [13]. Laser texturing can also reduce bacterial colonization on implants [15], while laser-induced periodic structures offer potential in encryption and anti-counterfeiting applications. Despite its indisputable interest, modelling self-organization of matter is complex because it involves electromagnetism, thermodynamics and fluid mechanics, whose interactions are still unknown. Consequently, there is currently no PDE that accurately models the whole

<figure><img src="image_1.png" /><figcaption>Fig. 1: Self-organization process triggered by ultrafast laser-matter interaction. On the right: emergence of a pattern after a series of double laser pulses.</figcaption></figure>

<footnote>This work has been partly funded by a public grant from the French National Research Agency (ANR) under the Investments for the Future Program (PIA), which has the reference EUR MANUTECH SLEIGHT-ANR-17-EURE-0026.</footnote>
physical phenomenon. Self-organization of matter faces a second issue: the impossibility to access abundant data. One of the major hurdles in acquiring large datasets stems from the inability of existing electronic devices for imaging on the fly the underlying dynamics going at the speed of light. The few data available is obtained through a complex setup involving a scanning electron microscope (SEM) to capture multiprocess images of an alloy sample after the laser-matter interaction has occurred.

Despite this unfavorable context, the ambitious goal of this paper is to better understand this complex dynamical system by leveraging the flexibility and efficiency of physics-informed machine learning. More specifically, we address the challenging prediction task of laser energy absorption of the surface. The energy distribution involved in this process is intricate, as it arises from a complex interplay between local nanoscale features. This phenomenon is governed by the unique characteristics of Maxwell's equations in inhomogeneous lossy mediums, requiring nano-scale physics modeling that comes with high simulation costs, and make any exploration impractical. A previous attempt to develop a surrogate for solving Maxwell's equations using physics-aware neural networks has been proposed in [7]. Known as MaxwellNet, this PINN-like model involves training a neural network by embedding the residual of Maxwell's equations into a physics-driven loss function. While our approach shares the focus on Maxwell's equations, it differs in significant aspects. MaxwellNet considers the propagation of a stationary light wave (we consider an ultra short laser pulse) through a transparent medium (we consider metals which involve light-matter interactions) in 2D (we consider 3D). Additionally, MaxwellNet predicts the electric field, while our focus is on the forecasting of the energy absorption resulting from the photoexcitation of free electrons. By directly predicting energy absorption, our aim is to reduce the computational time and gain quick insights into material responses across various electromagnetic environments, addressing multiple physics-related questions.

The contribution of this paper is three-fold: (i) We show that the way we pose the physical absorption problem allows us to overcome the data scarcity issue and train efficiently a simple surrogate model based on a U-Net architecture; (ii) Interestingly, our large experimental comparison study demonstrates that our model surpasses more sophisticated physics-informed neural architectures. We also show that our learned model allows a dramatic improvement compared to traditional Finite-Difference Time-Domain based solvers; (iii) An additional promising finding of our study comes from the correlation we exhibit between pattern formation and absorbed energy, prompting further inquiry into the still open question of the optimal arrangement of cavities/bumps favoring energy absorption.

The rest of the paper is organized as follows: Section II is devoted to the presentation of the preliminary back-

ground. Section III introduces our PiML model before performing a large experimental study in Section IV. We conclude in Section V and sketch promising lines of investigation.

## II. PRELIMINARY BACKGROUND

In this section, we recall the definition of Maxwell's equations as well as the Finite-Difference Time-Domain (FDTD) method used as numerical analysis technique. Even if this knowledge is not a prerequisite for understanding the rest of this article (eventually, we will simply solve an image-to-image regression task), reading the following sections A. and B. will give the reader a better grasp of (i) the underlying physics addressed in this study and (ii) the origin of our physics-driven simulation data that will be used for training the PiML models. Indeed, unlike Physics-informed Neural Networks (PINN [12]) that are trained while satisfying, through the minimization of PDE residuals, any laws of physics, our approach benefits from physical knowledge contained in simulation data. Like FNO [5] and U-NO [11], that are also briefly recalled in this section, our approach boils down to learning a mapping between two function spaces.

### A. Maxwell's equations and FDTD method

With advancements in computational power, the FDTD method, pioneered by [17], has experienced growing adoption for numerically solving Maxwell's equations. By incorporating the dielectric function $\epsilon_r$, the two curl operators in these equations can be reformulated as follows:

$$ \nabla \times \vec{E} = -\mu_0 \frac{\partial \vec{H}}{\partial t} \quad ; \quad \nabla \times \vec{H} = \epsilon_0 \epsilon_r \frac{\partial \vec{E}}{\partial t} + \sigma \vec{E}, $$

where $t$ is the time, $\vec{E}$ and $\vec{H}$ are the electric and the magnetic field, $\mu_0$ is the vacuum permeability, $\epsilon_0$ the vacuum permittivity, $\epsilon_r$ the relative permittivity, $\sigma$ the electrical conductivity.

In the FDTD method, a central differences scheme is employed within a Cartesian frame for both spatial curl operators and time derivatives. This approach facilitates the discretization of Maxwell's equations, enabling numerical solutions. The FDTD algorithm progresses through time using a leapfrog arrangement, wherein each component of $\vec{E}$ is updated and stored based on previous components of $\vec{H}$, while all components of $\vec{H}$ are updated from $\vec{E}$, and the cycle repeats. Additionally, to ensure numerical stability, the Courant condition must be satisfied. This condition dictates that $\Delta t < \frac{1}{c\sqrt{\frac{1}{(\Delta x)^2} + \frac{1}{(\Delta y)^2} + \frac{1}{(\Delta z)^2}}}$, where c is the speed of light in free space, $\Delta t$ the time step and $\Delta x$, $\Delta y$, $\Delta z$ are the spatial steps.

In summary, to numerically solve the equations using FDTD, it is essential to employ a central differences scheme for both spatial and temporal differentiations within a Cartesian coordinate system. This approach allows for the accurate representation of electromagnetic phenomena, including dispersive and lossy materials,
<figure><img src="image_2.png" /><figcaption>Fig. 2: Simulation domain. The numerical solver simulates the propagation of a laser pulse along the vertical axis and its 3D interaction with the surface roughness, and computes a measure of energy absorbed at every location and instant $A(x, y, z, t)$.</figcaption></figure>

while ensuring numerical stability through adherence to the Courant condition.

## B. Energy absorption

In this paper, we resort to the FDTD method to investigate the periodic energy deposition beneath a material's rough surface. The material features a roughness layer composed of cavities and bumps, as illustrated in Fig. 2. Positioned above the roughness layer is the source plane, situated within air, where a soft source is introduced into the computation grid. This soft source consists of an electric field component, either $\tilde{E}_x$ or $\tilde{E}_y$, representing an electromagnetic wave polarized along the $x$ or $y$ direction, respectively. This wave travels downward towards the material surface at normal incidence.

Throughout our simulations, we utilize linearly polarized light as the incident source, with the electric vector oriented either in the $x$ or $y$ direction. To minimize unwanted reflections at the boundaries, the computation domain is entirely enclosed by Perfectly Matched Layers (PML) [1], [2] in all six directions. PML are specialized boundary layers with evolving artificial conductivities. As depicted in Fig. 2, positioned around the boundary, these layers serve the crucial function of absorbing outgoing waves effectively, thereby minimizing reflections from the interface between the actual medium and the boundary medium. By keeping reflections from the actual boundary at acceptably low levels, PML play a vital role in simulating semi-infinitely extended media and preventing nonphysical reflections at the edges of the simulation grid.

The complex amplitude of the electric fields, denoted as $\tilde{E}$, is computed at intervals of half an optical cycle for every FDTD cell. This computation involves extracting the time-varying field value directly from the FDTD simulation. Subsequently, the energy absorption, denoted as $A$, is determined for each FDTD cell based on the complex amplitude of the electric fields, as follows:

$$A(x, y, z, t) = \frac{1}{2} c \epsilon_0 \frac{4\pi I m(\tilde{n})}{\lambda} Re(\tilde{n}) |\tilde{E}(x, y, z, t)|^2, \quad (1)$$

where $\lambda$ is the wavelength in free space, $\tilde{n}$ is the complex refractive index with $Re(\tilde{n})$ the refractive index (responsible for refraction) and $Im(\tilde{n})$ the extinction coefficient (describes absorption).

We implemented our own numerical solver for simulating the propagation of the laser pulse and computing the energy absorption from our experimental data. This solver will be used in the experiments described in Section IV.

## C. FNO and U-NO

Neural operators aim at learning a mapping between infinite dimensional spaces, including function spaces. They are known to be universal approximators of operators and have been shown to be very efficient for modeling dynamics in the form of PDEs. Two recent representatives of this family of models are Fourier Neural Operator (FNO) [5] and U-shaped Neural Operator (U-NO) [11]. As depicted in Fig. 3, FNO lifts (with L) the input (the surface topography in our case) to a higher dimension channel space and applies Fourier layers before projecting back (with P) and outputting the solution (the energy absorption in our case). Each Fourier block consists, on its upper part, in applying the Fourier transform followed by a linear transform on the lower Fourier modes before applying the inverse Fourier transform. To capture local information, a linear transform is applied in a lower branch (W). Finally, the composition of the two outputs is fed to a non-linear activation function.

On the other hand, U-NO also adopts the Fourier transform-based integration method, but differs from FNO by following a U-shaped architecture. It first performs an encoding by mapping the input to functions defined with smaller domains and then applies a decoding step to get a suited output function with skip connections from the encoder part. Like FNO, each element of the architecture takes the form of a Fourier layer.

Even though, at the first sight, both FNO and U-NO seem to be natural candidates for dealing with our task, we claim that they might face several cons for addressing the energy absorption prediction problem. First, the scarcity of data might prevent both methods from efficiently optimizing the large number of parameters. Second, as shown in the following section, the energy absorption prediction relies more on a image-to-image regression task than on the prediction of a dynamics, as targeted by FNO and U-NO. Finally, relying on a Fourier transform-based integration method, they might face issues for generating high frequencies.

## III. ENERGY ABSORPTION PREDICTION, SELF-ORGANIZATION DATASET, SIMULATIONS AND U-NET-BASED PML METHOD

In this section, we show how the prediction task of energy absorption can be dramatically simplified into a simple image-to-image regression task. This will have an important impact in terms of computational time and data
<figure><img src="image_3.png" /><figcaption>Fig. 3: Fourier Neural Operator (FNO) architecture.</figcaption></figure>

requirement for learning. Then, we describe the protocol we set-up for acquiring auto-organized surfaces. This involves fine polishing, irradiating and finally characterizing the surface by Scanning Electron Microscopy (SEM). This complex experimental setup is key for obtaining images that will be fed to the numerical solver in the experiments described in Section IV. Finally, we present the U-Net-based architecture used for learning an efficient surrogate solver allowing us to predict accurately and quickly the energy absorption.

### A. Towards an image-to-image regression task

Our goal is to predict the energy absorbed by the material when it is exposed to a laser pulse with a given polarization. For a given material, this mainly depends on the laser parameters (fluence, wavelength) but also on the topography that is typically represented as a binary 3D array, with 0 encoding for vacuum and 1 for the material. As introduced in Eq. 1, the absorption $A(x, y, z, t)$ is spatio-temporal. The sheer volume of data that this requires to manipulate and store is one reason why the numerical solver are costly to run.

To address this issue, we propose to redefine the prediction task so as to attain maximum training and prediction speed while keeping the physically meaningful information. This finally allows us to simplify the problem into an image-to-image regression task as described below. By cumulating the 3D topography along the propagation direction $z$, we obtain a 2D topography height map that can be viewed as a pseudo-3D array or 2.5D array since the sum managed to keep the third dimension information because the considered materials are not hollow. The output of our prediction task is also a 2D array, encoding the cumulated absorption at any $x, y$ position. More precisely and using an overloaded notation for the absorption, the task is to predict:

$$A(x, y) := \iint_{z,t} A(x, y, z, t) dz dt.$$

Integrating the temporal dimension $t$ is allowed by the fact that the laser pulses are ultrafast, more precisely at the femtosecond scale. Compared to other phenomena such as heating, melting and convection, the effect of a single laser pulse can be considered as instantaneous. Integrating along the light propagation direction $z$ is a compromise

to improve speed: we lose the information about the absorption at different depths but most of the complex effects happen in the first micrometers below the surface, with an exponential decrease with depth.

### B. Creation and acquisition of auto-organized surfaces

As our work fits in a global picture of understanding the formation of nanoscale self-organization patterns, we need to study the energy absorption on metal surfaces representative of the one used in self-organization physical experiments. These surfaces can be obtained through a process that requires a variety of skills and that we describe below. First, to achieve precise and high-quality laser processing of metal surfaces, surface preparation, involving an extremely fine polishing, is essential before laser irradiation. Next, the process involves irradiating the surface with a number of cross-polarized laser double pulses. Cross-polarized irradiation promotes self-organization by inducing isotropic energy deposition on the surface. It involves employing light beams with perpendicular polarization orientations to induce a specific effect. The latter is further enhanced by controlling the inter-pulse delay (in picoseconds), which governs the self-organization regime and the final pattern observed on the surface. To get a characterization of the surface topography after laser irradiations, SEM is employed to provide highly detailed images at significantly higher magnifications compared to conventional optical microscopes.

In this study, we consider 4 series of 51 SEM $800 \times 800$ images, some of them shown in Fig. 4. A polished sample is irradiated at 50 different locations with 1, 2, ..., 50 double pulses (a non-irradiated region is also considered, for a total of 51). The 4 series are obtained by varying the laser fluence and the delay between the two pulses (other parameters being fixed). We can observe that by changing these two parameters, we can obtain diverse patterns, including stripes, cavities or meshes. We can also note that after a certain number of double pulses, chaotic behaviors may appear. It is important to specify that these SEM images do not constitute our training set. To learn neural architectures from (input-output) pairs, we will need to perform simulation experiments from these images to obtain the corresponding (output) absorption (see Sec. IV).

### C. Neural architecture for predicting energy absorption

From the surfaces obtained (and imaged) through the process described in Sec. III-B, our approach is to learn to predict the 2D absorption as detailed in III-A. We will thus use the acquired SEM images to define our input topographies. Each topography will be fed to our numerical solver for Maxwell's equations (as explained in II-A). We will retain the 2D absorption $A(x, y)$ produced by the solver as our ground truth. Details about the training dataset generation are provided in Sec. IV-A. In terms of models, we suggest to learn a U-Net that is a go-to architecture for tasks that can be framed as image-to-image.
<figure><img src="image_4.png" /><figcaption>Fig. 4: SEM images of surfaces irradiated with cross-polarized double pulses, varying the number of double pulses between 0 and 50, for 4 combinations of laser fluence and delay (between the two laser pulses).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5: U-Net for predicting the absorption from the topography height map. All 3x3 convolutions are followed by a ReLU activation and a batch normalization layer.</figcaption></figure>

U-Net is a convolutional network originally designed for semantic segmentation tasks, particularly in biomedical image analysis [14]. This architecture allows the network to capture both local and global information efficiently, making it well-suited for tasks where precise localization is essential. However, its design can be repurposed for regression tasks to predict continuous values instead of discrete labels. Thanks to the simplification of the energy prediction task (as described in Sec. III-A), we can now leverage this property of U-Net. The regression neural architecture, depicted in Fig. 5, comprises two key components: the contracting path, which extracts features and reduces spatial dimensions, and the expansive path, which restores spatial resolution for precise localization. The former includes double convolution blocks for feature extraction and downsampling blocks for spatial reduction. The latter involves upsampling blocks to restore resolution and double convolution blocks for refining features. All 3x3 convolutions are followed by a ReLU activation and a batch normalization layer. The final output layer, a single convolutional layer with a kernel size of 1, aggregates features from the expansive path to produce the regression output which is fed to a MSE loss function (see Eq. 2).

## IV. EXPERIMENTS

In this section, we first present the training dataset specifically generated for this study. Then, we detail the learning setup and the PiML methods compared, before analyzing the results.

<figure><img src="image_6.png" /><figcaption>Fig. 6: Illustration of examples of the dataset composed of pairs of (topography, absorption) images. We use rotated images so that all examples virtually have the same polarization.</figcaption></figure>

### A. Simulation data generation from SEM images

To utilize SEM images and construct the dataset using our numerical solver employing the FDTD method, we made several assumptions. Given that our solver operates in a 3D spatial domain to compute solutions of Maxwell's equations, we transformed our scanned 2D images into 3D representations. This transformation involved adding depth under the assumption that brighter pixels in the image are closer to the SEM sensor than darker pixels. To ensure accuracy, we supplemented the surface with additional material underneath, accounting for the PML to dissipate sufficient energy. In our simulations, we begin by establishing an orthonormal reference frame within the provided space. From there, we simulate the propagation of a wave traveling downward towards the material surface at normal incidence along the *z*-axis. Throughout these simulations, we employ linearly polarized light as the incident source, with the electric vector aligned either in the *x* or *y* direction to derive $A(x, y) = \iint_{z,t} A(x, y, z, t) dz dt$. For a given surface, the numerical solver can actually generate 2 different values of $A(x, y)$. As illustrated in Fig. 6, when nanostructures are present, the pattern of absorption strongly depends on electromagnetic wave polarization. Let us denote $S(x, y, z)$ a surface before laser irradiation and, with an additional notation overloading $A_0(x, y)$ (resp. $A_{90}(x, y)$) the energy absorption respectively for electric vector aligned in the *x* (resp. *y*) direction. In order to enforce the equivalence assumption, which is more data-efficient, we rotated the initial surface to generate $A_0(x, y)$ and $A_{90}(x, y)$ from $S(x, y, z)$, as illustrated in Fig. 6.

We conducted numerical simulations involving the four sets of 51 scanned surfaces already depicted in Fig. 4, each characterized by a refractive index of $n = 2.15$ and an extinction coefficient of $k = 4.3$. To reduce SEM noise, we preprocess these images using a Gaussian blur filter with $\sigma = 4$. We use our numerical solver to simulate the irradiation by a femtosecond laser with a wavelength of $\lambda = 1.03\mu m$, polarized along either the *x* or *y* axis. To work at physically relevant scale and fit into the hardware constraints (GPU memory), the simulation spatial resolution is chosen as $dx = dy = 1.67nm$ and $dz = 1nm$. The spatial domain covers roughly
<table><caption>TABLE I: Comparison of the 5 PiML methods on the 4 test series in terms of errors (mean and standard deviation σ) as well as average training time and number of parameters.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">Test (0.18_8)</th><th colspan="2">Test (0.21_2)</th><th colspan="2">Test (0.25_36)</th><th colspan="2">Test (0.26_28)</th><th rowspan="2">Av. training time (s)</th><th rowspan="2">Nb of parameters</th></tr><tr><th>mean</th><th>(σ)</th><th>mean</th><th>(σ)</th><th>mean</th><th>(σ)</th><th>mean</th><th>(σ)</th></tr></thead><tbody><tr><td>U-Net</td><td>8.6</td><td>(9.9)</td><td>6.0</td><td>(6.5)</td><td>5.1</td><td>(4.5)</td><td>12.0</td><td>(14.6)</td><td>1740</td><td>~2e+6</td></tr><tr><td>FNO</td><td>20.2</td><td>(22.3)</td><td>13.7</td><td>(16.0)</td><td>10.0</td><td>(11.3)</td><td>19.2</td><td>(20.3)</td><td>5127</td><td>~1e+8</td></tr><tr><td>FNO conv3x3</td><td>19.3</td><td>(23.8)</td><td>14.4</td><td>(15.0)</td><td>9.6</td><td>(10.8)</td><td>18.9</td><td>(20.9)</td><td>5247</td><td>~1e+8</td></tr><tr><td>FNO+U-Net</td><td>16.7</td><td>(20.7)</td><td>11.2</td><td>(13.9)</td><td>8.9</td><td>(9.3)</td><td>17.7</td><td>(20.0)</td><td>13317</td><td>~1e+8</td></tr><tr><td>U-NO</td><td>25.2</td><td>32.4</td><td>15.7</td><td>18.5</td><td>10.9</td><td>11.9</td><td>30.8</td><td>33.2</td><td>4273</td><td>~1e+8</td></tr></tbody></table>

1.5μm, with the material occupying half of the z direction, including 20nm of roughness. The PML size is set to 50 (cells). With two polarizations, we generated from the solver pairs of topography and absorbed energy, resulting in a total of 408 = 4 × 51 × 2 pairs. To expand our dataset and learn on smaller patches, we extract 9 crops per image using a sliding window scheme to create new pairs of size 320×320 pixels, leading to a final count of 3672 = 408×9 training pairs (S, A) with an overlapping. It is important to note that these simulations are purely computational and do not involve physical experimentation. Furthermore, we normalized the 2D energy absorption by the energy absorbed by a flat surface. This normalization is essential because the numerical solver uses a source with a non-uniform fluence and was not calibrated to provide meaningful output in terms of raw values. By normalizing, we can correct this issue, allowing this quantity to be interpreted as an increase or decrease in energy absorption compared to a flat surface.

## B. Experimental setup and PiML methods compared

In this study, we compare 5 different PiML methods: (i) **U-Net** as described in Sec. III-C with the number of (down/up)-sampling blocks ∈ {3, 4, 5, 6} as hyperparameter; (ii) **FNO** [5] with the following hyperparameters: number of Fourier layers ∈ [3, 4, 5], the size of the lift L ∈ [16, 32, 64], the number of Fourier modes ∈ [128, 192, 256]; (iii) **FNO conv3x3**, a variant of FNO with a convolution W of size 3x3; (iv) **FNO+U-Net**, a variant exploiting the composition of the outputs of FNO and U-Net; (v) **U-NO** [11] with the optimal hyperparameters obtained with FNO and `uno_out_channel=[16,32,32,16]`, `uno_scalings=[[1.0,1.0], [0.5,0.5], [0.5,0.5], [4.0,4.0]]`. For each method, we repeat 4 times (once for each series playing the role of a test set) the following leave-one-series-out cross-validation procedure: 2 series of SEM images are used to train the models with all hyperparameters configurations, the third series being employed to evaluate the validation accuracy; the process is repeated 3 times over the 3 folds and the average validation accuracy is used to decide on the best hyperparameters; Finally, the latter are used to learn a model from the 3 series and evaluated on the test series.

For each method, the training process consists in learning a model $h_{\theta}^*$ that outputs a 320 × 320 image of energy absorption given a topography S. $h_{\theta}^*$ is parameterized by $\theta$

(the parameters of the considered architecture) and learned by solving the following problem:

$$h_{\theta}^* = \arg \min_{h_{\theta}} \frac{1}{N \times 320^2} \sum_{i=1}^{N} \|h_{\theta}(S_i) - A_i\|_{\mathcal{F}}^2, \quad (2)$$

where N is the total number of pairs $(S_i, A_i)$ representing the pairs of (topography,absorption) images and $||.||_{\mathcal{F}}$ is the Frobenius norm. In our leave-one-series-out procedure, $N = 2 \times 51 \times 2 \times 9 = 1836$ for each run. Note that we solve Problem 2 with a weight decay equal to $1e-7$. All the networks have been trained for 100 epochs using an Adam optimizer.

## C. Analysis of the results

The results are reported in Table I where the mean (absolute error) and the standard deviation $\sigma$ are computed over all the pixels of the test images. In order to interpret the magnitude of these numbers, it is worth noting that the range of values for the energy absorption depends on the topography of the pattern. As illustrated in Fig. 7 (second column), it can vary from 0 to about 1000. Taking into account this information, we can see that all the PiML methods perform pretty well for predicting the energy absorption. But the most striking remark we can make from this table is that our simple image-to-image U-Net architecture outperforms the other competitors while being much cheaper in terms of computational time and number of parameters. This finding confirms our intuition that the other PiML approaches struggle to capture the dynamics with such scarce data. We can also notice that the model behaves well on the four series even if predicting the last one (0.26_28) seems more challenging. In this low data regime, this can be explained by the fact that this fourth series is mainly composed of cavities that are underrepresented in the other series. Finally, we can note that the 5 models work very well on the third series. This can be explained by the fact that this latter does not exhibit any chaotic behavior with a late appearance of the first patterns, leading to a much simpler prediction task.

We illustrate in Fig. 7 some representative prediction results from the four series. The visual inspection confirms the ability of U-Net to generalize to unfamiliar patterns. The second (Ground truth A) and third columns (Predicted $h_{\theta}(S)$) illustrate the capacity of U-Net to outputs very similar images. The relative difference $(h_{\theta}(S) - A)/A$ is
<figure><img src="image_7.png" /><figcaption>Fig. 7: Illustration of prediction results. From left to right: Topography S, ground truth energy absorption A, predicted energy absorption $h_θ(S)$ by U-Net, relative error of U-Net prediction $(h_θ(S) - A)/A$, relative error histogram, relative error of FNO, FNO conv3x3, FNO+U-Net and U-NO.</figcaption></figure>

reported in the fourth column, where the lighter the better. To emphasize the good prediction ability of U-Net, we report on the right hand side of Fig. 7 the relative error for the four competitors, exhibiting much more pronounced differences, yet still acceptable. The success of U-Net despite scarce data can be attributed to the dimensionality reduction we imposed. By transforming the complex 3D topography into a 2D height map that retains essential height information, the model can more efficiently learn and identify critical features that influence energy absorption. This simplification is particularly advantageous because a 3D-to-2D task typically requires more data to capture features accurately.

### D. Speedup and interactive exploration

The numerical solver used is implemented on GPU. Our simulations have been sized based on the GPU memory available to run the simulations (24GB). The average runtime for a simulation on 1 SEM image is just above 30 minutes (exactly 1857sec) which corresponds to 8.5 GPU-days for the 408 images.

Training the network for 100 epochs takes less than 30 minutes. The time to process a full 800 × 800 image with our U-Net (the network being fully convolutional, it is not limited in the image size it can take as input), on the same machine and averaged over the whole set and with a batch size of 30, is 12.5ms per image. Compared to the 30 minutes, this corresponds to a speedup over $10^6$ (exactly 144000×). With a batch size of 1, simulating the case where we want to process images as soon as possible, the latency goes up to 97ms in average. The CPU version takes 1.3 second per full-size image. Given the speedup of the model, we created an offline web application that allows the user to explore interactively the absorption predicted by the model, working with the nominal input size for our U-Net, namely 320×320 images<fnref n="1" />.

### E. Exploring physics questions using the fast proxy

We consider a couple of physical interrogations around energy absorption and analyze whether conclusions drawn with our U-Net would be coherent with the ones drawn from the numerical solver. Both questions are concerned about the total absorption of a topography, which is a scalar value obtained by integrating the absorption images spatially and temporally as follows:

$$A := \iiint_{x,y,z,t} A(x, y, z, t) \, dx \, dy \, dz \, dt.$$

The initial inquiry pertains to the overall absorption's evolution as the number of double pulses increases, taking into account both pulses with polarizations at 0 and 90 degrees, denoted as $A_0 + A_{90}$. To improve interpretability, we normalize the total absorption by its initial measurement (at $N = 0$ impulsion), as the numerical solver was not calibrated to produce meaningful output in its original scale. Predicted and ground truth overall absorption are plotted in Fig. 8 showing the capacity of the model to predict the absorption tendency, yet having a little systematic bias, i.e. under-estimating the absorption.

Initially, during the first few double pulses, the actual energy absorption remains relatively consistent and low. However, after numerous pulses, there is a significant increase, eventually reaching a plateau. Sometimes, as the material undergoes multiple irradiations, the energy absorption decreases due to deterioration. Very interestingly, the absorption keeps increasing for the third series. Observing this series from Fig. 4, we can conclude that this definitely makes sense as the dynamics has not yet reached a chaotic behavior. The predicted energy absorption follows a similar trend to the actual absorption, enabling us to draw the same conclusion.

A second question is about the relative absorption between both polarizations. More precisely, we are interested in determining whether the self-organized structures exhibit differential absorption between the two polariza-

<footnote n="1">¹The code of the web application as well as that of the experiments and the datasets will be made available upon acceptance of the paper.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 8: Overall absorption (sum of the two pulses) as a function of the number of cross polarized pulses within the series. The absorption tends to increase with the self-organisation. The prediction follows the ground truth.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 9: Absorption difference between the two polarizations presented as a ratio (1 means no difference). The first two series show no directional preference while the others get specialized to an orientation after self-organization starts. The prediction follows the ground truth.</figcaption></figure>

tions or if they are direction-independent. Therefore, we computed and plotted the ratio between $A_0$ and $A_{90}$. Fig. 9 shows the predicted and ground truth ratios, as a function of the progressive auto-organization of the surfaces.

Here again, the model prediction follows the same tendency as the ground truth from the simulator and would lead to the same overall conclusions, being it for the two series that exhibit orientation specialization and the other two that show no specificity (values oscillating around 1).

V. CONCLUSION AND PERSPECTIVES

In this paper, we efficiently deal with the prediction of energy absorption of irradiated surfaces by simplifying the task into an image-to-image regression problem. We show that in the low data regime characterizing the application at hand, learning a U-Net architecture leads to better results than more sophisticated PiML models. The experimental results show a good accuracy and exhibit a promising correlation between the absorbed energy and the self-organization of patterns. Our main perspective involves using this very fast and differential proxy of the absorption as a building block for end-to-end learning of the (non-observable) dynamics that happens between and after the two pulses. The control of the self-organization process upon multi-pulse excitation is inherently complex, as a dissipative structure adheres to laws governing energy and entropy flux optimization that are not fully understood. In these delicate conditions, mastering laser parameters to achieve a uniform and reproducible pattern is imperative for precise surface engineering. Predicting and understanding absorption mechanisms is thus the crucial first step in guiding the structuring process during irradiation. Access to electromagnetic absorption for any type of surface

opens up the possibility of considering a much broader range of conditions, thereby enabling more efficient design of metasurfaces essential for the development of optical devices, sensors, and energy-harvesting systems. Potential applications of the model include its ability in inverse design tasks, such as identifying structures optimized for maximal or minimal absorption, or with specific spatial distributions tailored to particular requirements. Additionally, restricting predictions solely to 2D absorption may limit exploration of complex phenomena. Thus, extending predictions to 3D absorption, either at the simulator's scale or, for efficiency, at a coarser resolution (e.g., predicting three absorption values per *x*, *y* location: surface, shallow, and deep) may be necessary to capture the full scope of the phenomenon under study.

REFERENCES

[1] Bérenger, J.P.: A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics 114(2), 185-200 (1994)

[2] Berenger, J.P.: Three-dimensional perfectly matched layer for the absorption of electromagnetic waves. J. Comput. Phys. 127(2), 363-379 (sep 1996)

[3] Hao, Z., Liu, S., Zhang, Y., Ying, C., Feng, Y., Su, H., Zhu, J.: Physics-informed machine learning: A survey on problems, methods and applications (2023)

[4] Karniadakis, G.E., Kevrekidis, I.G., Lu, L., Perdikaris, P., Wang, S., Yang, L.: Physics-informed machine learning. Nature Reviews Physics 3(6) (5 2021)

[5] Li, Z., Kovachki, N.B., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A.M., Anandkumar, A.: Fourier neural operator for parametric partial differential equations (2021)

[6] Li, Z., Zheng, H., Kovachki, N., Jin, D., Chen, H., Liu, B., Azizzadenesheli, K., Anandkumar, A.: Physics-informed neural operator for learning partial differential equations (2023)

[7] Lim, J., Psaltis, D.: MaxwellNet: Physics-driven deep neural network training based on Maxwell's equations. APL Photonics 7(1), 011301 (01 2022)
[8] Long, G., Lu, W., Zha, L., Zhang, H.: PHC-GAN: physical constraint generative adversarial network for single image dehazing. In: 32nd IEEE ICTAI conference. pp. 545–549. IEEE (2020)

[9] Long, Z., Lu, Y., Ma, X., Dong, B.: Pde-net: Learning pdes from data (2018)

[10] Ozan, S., Bilgin, A., Şefika Kasman: Laser textured ti-6al-7nb alloy for biomedical applications: An investigation of texturing parameters on surface properties

[11] Rahman, M.A., Ross, Z.E., Azizzadenesheli, K.: U-NO: U-shaped neural operators. Transactions on Machine Learning Research (2023), https://openreview.net/forum?id=j3oQF9coJd

[12] Raissi, M., Perdikaris, P., Karniadakis, G.E.: Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics (2019)

[13] Ranella, A., Barberoglou, M., Bakogianni, S., Fotakis, C., Stratakis, E.: Tuning cell adhesion by controlling the roughness and wettability of 3d micro/nano silicon structures. Acta Biomaterialia 6(7), 2711–2720 (2010)

[14] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation (2015)

[15] Sarbada, S., Shin, Y.C.: Superhydrophobic contoured surfaces created on metal and polymer using a femtosecond laser. Applied Surface Science 405 (2017)

[16] Vorobeyev, A.Y., Guo, C.: Colorizing metals with femtosecond laser pulses (2008)

[17] Yee, K.: Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media. IEEE Transactions on Antennas and Propagation