On the use of a digital twin to enhance femtosecond laser inscription of arbitrary phase patterns

Thanks to the non-linear nature of laser-matter interaction, the use of femtosecond lasers offers a versatile method for encoding information and modifying transparent materials in their volumes and this, with sub-micron resolution. The underlying physical process is a succession of intricate and complex nonlinear phenomena that are sensitive to multiple and multidimensional parameters, such as beam intensity distribution, exposure dose homogeneity and pulse-overlapping sequences as well as propagating wavefront angular orientations and temporal distortions. As a consequence to this inherent and often overwhelming complexity, obtaining a repeatable and accurate result relies strongly on time-consuming machine-specific calibration and experiment-specific fine-tuning attempts until the desired result is reached. Here, we present a digital twin of the processed specimen that not only accurately predicts the exposure outcome in terms of introduced retardance, but also offers a pathway for designing feedforward schemes that compensate for known inaccuracies. We demonstrate the merit of this approach through illustrative examples of arbitrary phase patterns, forming waveplates and images, based on refractive index modulation induced during laser exposure.


Introduction
Tightly focused ultrashort laser pulses allow for non-ablative structural modifications in the bulk of dielectric materials. These modifications take different forms depending on processing parameters (e.g. pulse energy and duration, scanning speed, repetition rate) and processed materials [1]. In dielectrics, several technologies exploit these effects, most notably direct waveguides inscription [2] or inscription of channels of increased etching rate, allowing for versatile three-dimensional micro-fabrication [3][4][5].
The full process leading to structural modifications is a succession of several phenomena, leading to various outcomes depending on complex relationships between process parameters [6]. These processes are highly sensitive to subtle differences in exposure conditions [7] and interaction between neighbouring exposed regions due to stress build-up in the material [8,9]. A dramatic manifestation of the inherently complex and rich phenomena underlying non-ablative femtosecond laser exposure is found in the anisotropic behaviour of the exposure process observed in otherwise isotropic materials [10][11][12]. There, radically different patterns are produced by simply reversing the laser writing direction, despite identical pulse energy and writing conditions. Consequently, apparently benign changes in the laser source or through the beam path may have a direct impact on the process reliability, as the pulses' characteristics might change dramatically. Furthermore, several sources of disturbances and inaccuracies are inherently present in laser-based manufacturing platforms, may they be erratic or repeatable. As the laser is scanned over the specimen (or likewise the specimen moved under the laser beam) to define a trajectory, for a set pulse repetition rate, the energy exposure dose-a key parameter defining the process outcome, varies during acceleration and deceleration phases due to intrinsic dynamical limitations. Other significant sources of static disturbances, particularly difficult to correct for, are found in geometrical (ellipticity) and spatio-temporal characteristics (spatial chirp, or pulse-front tilt [10,11,13]) of pulses themselves. Although these issues may be mitigated using either extensive calibration, adjustment of the optical systems, or by adding adaptive optics, this comes at the cost of highly complex laser platforms with limited throughput performances and robustness. These technical constraints of optical and mechanical origins make the process path-dependant. Consequently, two trajectories leading to the same pattern may not yield a priori the same phase map as will be shown in the sequel.
Here, we propose a general framework relying on a digital twin to predict and enhance the inscription of phase maps embedded in dielectrics, for arbitrary scanning strategies and this, without adding further complexity to the laser inscription setup. As a twin of the real experiment, the digital twin offers the possibility to test the outcome of arbitrary exposure patterns, without effectively performing any physical experiments and can be used to test and design feedforward algorithms for direct-inscription of precise arbitrary phase maps with high throughput and accuracy.

General principle
A digital twin is a virtual replica of a physical system that mimics its behaviour in a virtual world. In the context of laser inscription, it is a virtual model of the material response to laser exposure that simplifies optimisation of process parameters, alleviating the need for extensive and pattern-specific calibration and that offers a virtual framework for testing exposure patterns. The ideal digital twin embeds all the physics of a system up to a point where one cannot truly distinguish the differences between virtual and experimental outcomes. Reaching this level of replication is an ultimate endeavour. In this illustration, the digital twin outputs a single observable quantity: a map of laser-induced retardance, used as a process metric. It inputs exposure conditions (i.e. pulse, processing parameters and beam path) and outputs a two-dimensional phase map.
In glass systems, femtosecond laser illumination leads to localised changes of optical retardance that have different physical origins depending on the material considered. It may originate from density changes, from form-birefringence due to sub-optical-wavelength density modulation (as first observed in fused silica [14]), as well as by localised composition changes or crystallisation effects. Due to the small size of the laser modified zone (not exceeding a few microns at most along the optical axis), laser-induced retardance changes remain small, and effectively measuring them requires complex optical setups, difficult to implement for in situ observation. Here, we use recent progresses in digital holographic microscopy (DHM) to probe changes in optical retardance directly on the processing platform.
The digital twin takes into account imperfections related to the specimen positioning under the laser beam using prerecorded effective positions for a given pattern. This ensures that the local exposure dose is accurately estimated. To simulate a laser-exposure experiment, the energy modulation signal is calculated according to the desired pattern output and synchronised with the trajectory. Figure 1 illustrates this working principle for a simple trajectory consisting of concentric squares.

Modelling of the effect of exposure dose and beam anisotropy
The predictive model that captures both beam imperfections and the effect of the exposure dose on the material is inspired from the general relationship devised by Chichkov et al [15] for ablation, itself adapted from the two-temperature model [16], that we adapt to the case of non-ablative dose accumulation.
In polar coordinates, the retardance induced by a pulse set on (0, 0) as a function of the maximum fluence F, the apparent fluence threshold F th , and a polar beam energy distribution Ξ(ρ, θ) writes: ) . (1) F th is a threshold fluence, defined here as the energy for which one pulse will induce a measurable refractive index modification of the material. r 0 is a normalised retardance constant, in this study set to unity, as we do not claim to model the exact value of retardance, but its variations. In what follows, we consider 'pulse packets' rather than individual ones to reduce computation time given the high repetition rates used, ν rep = 100 kHz. As the motion controller retrieves positions at a steady but lower frequency in the range of 100-500 Hz, these points are directly used to define these pulse packets in time and space. This assumption is reasonable considering the gradual changes of index upon accumulation of pulses. Given the speed v scan at which the laser beam moves in the specimen, consecutive pulses overlap significantly (>99%). In practice, F th is set to unity as the fluence threshold is calibrated before inscribing. As the nth packet influences the following ones, the influence of previously written structures needs to be considered. We take into account a dependency of the value of F th with the local retardance before absorption R n − 1 as follows: with α being an empirical factor. Ξ(ρ, θ) is a polar distribution of the modification profile. This element portrays the beam energy profile at the focus as well as possible sources of anisotropicity, such as beam ellipticity or pulse-front tilt. In a first approximation, we encapsulate all these effects in a single modified Gaussian distribution according to: with e being the ellipticity of the distribution, θ 0 its absolute orientation, and w 0 a length equal to the isotropic standard deviation in the e = 0 case. As this distribution embeds multiple phenomena (and not only a possible beam-intensity elliptical distortion), both parameters, e and θ 0 are devised empirically. This distribution is assumed to be valid in the case of a Gaussian beam profile. Finally, we rewrite F/F th as (1 + β) to simplify further. By developing all elements in the expression of R n (ρ, θ), we finally obtain the expression on equation (4).

Discrete pattern description
To compute the digital twin output, we define a 2D array that samples the inscription plane in a polar coordinate system. This array is named R n [z, z n ], with n being the packet's ordinal number and z the discrete real position of the array in a complex form, and z n the central position of the packet. To further reduce computing time, each iteration is calculated locally, i.e. in the vicinity of the central position of the packet (square of 2w 0 of side length). Here, we encode the pattern by locally modulating the energy that will yield to different retardance level. To account for it, the normalised fluence term F/F th or 1 + β is made variable as a function of the packet ordinal n such as 1 + β n . As initial condition, we consider a homogeneous and entirely pristine substrate, therefore ∀z, R 0 (z) = 0. Finally, equation (4) becomes: Model parameters are optimised by comparing the non-corrected model to the process output, similarly to the case shown below in figure 5. The nonlinear, bound-constrained limited-memory  (2) and (3) is evident as this beam path exhibits a significant path dependence.
Broyden-Fletcher-Goldfarb-Shanno optimisation algorithm (L-BFGS-B) [17] is used to fit the different parameters to minimise the mean squared error between the two maps. This optimisation step is performed for any major change in experimental conditions (e.g. beam profile, pulse duration).

Determination of a modulation signal to compensate for scanning static inaccuracies and beam asymmetries
Let us now consider the case of adding a feedforward algorithm to compensate for predicted deviations from the actual result and the expected one. To implement such an algorithm, one could run several experiments and optimise the outcome based on the measured results. Here, the digital twin provides a faster path forward as it provides a digital and accurate replica (as will be further demonstrated in the following section) of what is to be expected should the physical exposure be done. To illustrate of possible feedforward algorithms, here we consider a first order polynomial dependence on the measured difference between the expected pattern and its predicted version.
More complex algorithms, including non-linear dependencies could be considered, following a similar methodology.
The values of β n are calculated from the expected retardance level as defined by the pattern, expressed according to the laser path, R des [s n ] (s n being the position of the packet on the beam path) and the local deviation of the model output from a base value corresponding to the case without imperfections ∆R N [s n ]. The latter corrects for the aforementioned imperfections in the inscription process. The following equation (6) expresses the modulation signal β n .
with γ 0 , and γ 1 being parameters optimised by comparing corrected model outputs with the inputs. Here, we use the BFGS algorithm [17] to perform this task. This optimisation step is needed for any change in model parameters. As an illustration of outputs given by the model, we show one specific path strategy (concentric squares) and the corresponding output on figure 2. The difference in the case of the uncorrected outcome is significant.

General process
In this experiment, the specimen is moved under the laser beam using translation stages (PI Micos GmbH). The use and implementation of the digital twin is illustrated in figure 3. The stages' repetitive tracking errors over the scanned trajectory are taken into account by acquiring the effective stage positions, as recorded by the high accuracy linear encoders, during a calibration run with the laser off at a sampling rate of 100-1000 s −1 . The sampling rate is adjusted to avoid buffer overflow, with a maximum of 120 samples per individual move. These inputs are uploaded into the digital model described previously, which define a Then, a greyscale image defines the pattern, and therefore the desired output. Using both, a pulse energy modulation signal is computed for the selected scan path and applied using an acousto-optic modulator (AOM). The digital twin predicts the expected retardance map as the outcome of the process. The modulated signal serving as input to the digital twin can be a one-to-one relation with the desired pattern or a more complex modulated signal that integrates feedforward compensation. As an illustration, the illustrative pattern is a prehistoric painting of horses-one of the oldest human representation from prehistoric times dated back from 30 000 BCE-from the Chauvet cave at Vallon Pont-d'Arc (France).

Figure 4.
Processing setup with a transmission DHM (from Lyncée Tec) sharing the same beam path as the laser beam used for inscribing the pattern. The femtosecond laser beam is provided by an Yb-based regenerative amplifier, modulated using a free-space acousto-optic modulator (AOM). The beam size is adjusted using a beam expander (BE) to match the pupil of the objective. The initial power level is set using a rotating half-wave plate (λ/2) and a polarised beam-splitter (PBS), and the direction of the linear polarisation state using another half-wave plate. A fibre pigtailed 794.5 nm-CW diode laser source is used for illumination. The beam is split between an object beam, slightly defocused to move its effective focal spot relative to the injection objective, and a reference beam. Both are merged to generate a hologram on the camera. The PC acquires and processes the image to obtain the phase map. It further drives the motion controller as well as the AOM through a digital-to-analogue converter (DAC) commanding an RF driver. The modulation signal is derived from the digital twin. modulation signal, according to the desired output given in greyscale bitmap format. The modulation signal is then uploaded into a digital-analogue converter (DAC), which drives an RF-driver that commands the acousto-optic modulation of the laser pulse intensity. To ensure proper synchronisation of the translation of the stages with the modulation signal, the DAC and the controller are synchronised using the controller as the master. The model parameters are tuned using in situ phase-contrast microscopy, using a DHM from Lyncée Tec, integrated on the processing platform as shown in figure 4. The model relies on a limited set of parameters, and given the relatively fast execution of the algorithm (5-20 s on a 2.4 GHz Intel Core i5 processor), manual adjustment of the model parameters is performed in a few minutes.

Processing platform
The femtosecond laser source used for processing is an Yb:fibre mode-locked oscillator, coupled at 100 kHz into an Yb-based regenerative amplifier (Amplitude Systèmes SA). This system delivers 150 fs FWHM pulses at 1030 nm (GRENOUILLE measurement, Swamp Optics). An acousto-optic modulator diffracts the output beam, allowing fast amplitude modulation of the pulse train. The beam is expanded using a 2-5× adjustable beam expander to fit the pupil diameter of the objective microscope, in order to reach the maximum numerical aperture (NA). A half-wave plate controls the linear polarisation state of the beam. The beam is then focused using a 20× magnification, 0.40 NA, high-power objective lens (Thorlabs Inc. MicroSpot®line, the average spatial sampling period). In what follows, we used fused silica square substrates (Corning 7980 0F) as test specimens.

Phase-contrast imaging system using an on-axis DHM
In situ monitoring of the phase map for a region of interest (ROI) in the vicinity of the laser inscribing focal point is performed using an on-axis DHM (Lyncée Tec, Switzerland). This technology allows us to bypass the need for complex optical elements to retrieve the effective retardance. The source used is a fibre pigtailed CW diode functioning at 794.5 nm. The fibre is split between an object beam path, that will perform bright field imaging in a conventional manner, and a reference beam path used to generate the hologram [18,19]. The object beam is collimated from the fibre-end and is merged with the femtosecond processing beam path using a dichroic mirror. Before the latter, a 250 mm plano-convex lens slightly defocuses the beam before the processing objective lens, set in order to reduce its effective NA and displace the actual focal point. The object beam therefore irradiates an area larger than the ROI for imaging purposes, while the processing beam is diffraction limited. A differential interference contrast objective lens of 100× magnification images the ROI. The reference beam is then collimated and merged with the object beam after the imaging objective lens to generate the hologram acquired on a CMOS camera (typically used at 10 fps). The reconstruction algorithm is running in quasi-real time on a conventional personal computer (PC). It decomposes the hologram into an amplitude and a phase map, the latter being the quantitative phase-contrast image used as a metric for the process outcome and for comparing with the predicted phase map that the digital twin outputs.

Experimental demonstration and discussion
First, we illustrate a case where the pattern orientation is changed and how the digital twin is able to predict the correct outcome. This example underlines the effect of beam anisotropy on the outcome of a laser exposure. Second, we show how the digital twin can be used to correct for inhomogeneity using the feedforward algorithm described earlier.

Effect of pattern rotation
We consider a simple square pattern (40 µm × 40 µm, shown in figure 5) of homogeneous refractive index as illustrative example. As the pattern is one order of magnitude larger than the beam size itself, scanning-at a given sampling level-is required to achieve the desired result. Clearly, there exist multiple possible trajectories that theoretically fulfil the same outcome. The choice for one or another will rely on technological implementation constraints, such as-but not only-the scanning system dynamics.
Here, for illustrative purpose, we choose a scanning trajectory consisting of a concentric pattern with vertical and horizontal edges, written clockwise. A pitch of 0.5 µm is chosen, which corresponds to about a third of the linear spot size. Such pattern is non-optimal as it leads to severe acceleration and speed constraints on the moving stages as they near the centre of the pattern. This choice was made to test the accuracy of our digital twin and to challenge its robustness for capturing anisotropic effects. The same pattern is written for different orientations from a common axis. Although a homogeneous pattern would have been theoretically expected, four sectors are revealed depending on the pattern orientation from the experimental phase maps. The strong dependency on the pattern orientation is remarkable and is a direct consequence of the preferred orientation of the incoming beam. Note that this difference is relatively unnoticeable at a 45 • rotation angle. Although the laser was optimised to correct for ellipticity in the intensity distribution and despite the material is itself isotropic (fused silica), the revealed preferred-orientation relates to additional inherent beam anisotropicities that could not easily be corrected for, such as pulse-front tilt.
In our digital twin, we consider an ellipsoidal distribution of the laser-induced retardance to capture such anisotropic effects, with the long axis oriented at a given angle of 0 • . As can be seen in figure 5 that compares the digital twin predictions with the output of real experiments, the digital twin gives a realistic rendering of the actual experimental observations.   France). Note that the DHM is used at the limit of its effective resolution and consequently, the improvement appears somewhat limited. The true resolution of the laser-written patterns is effectively higher. This point is addressed in section 4.4.

Inscription of arbitrary patterns
To validate the feedforward optimisation scheme, we consider a set of arbitrary patterns, as shown in figure 6(b). We selected the worst-case scenario corresponding to the maximum anisotropy observed in figure 5, i.e. a rotation angle of 0 • with respect to the optical axis and for a non-optimal trajectory that involves acceleration and deceleration sequences. The results are shown in figures 6(c) and (d). Although the conditions are non-optimal and quite detrimental, one can see that even in such 'harsh' conditions, the algorithm is able to correct for the main anisotropic effects. The pattern edges still appear significantly in the second case where the material should be pristine (e.g. third and fourth case): the high density of points there, although at a low fluence threshold, is enough to generate some excess retardance. On the other hand, we observe a slight overshoot of the modulation where the substrate should be inscribed (first and second case): there the model predicts a stronger retardance.
Our model does not consider the effect of stress inherently present with laser-written patterns [9]. Taking it into account could further improve the quality of the results, in particular when pattern discontinuities are present and high density of lines are used. In the case illustrated here, this effect remains small as the laser-pulse duration and pulse energy used here induce localized densification in the material. In this particular regime, the laser-induced stress is small and can be made minimal by tuning the net exposure dose [20]. Nevertheless, enriching the model in that regard is an interesting outlook. This could be achieved by including a separate map describing the variations of density in the vicinity of a laser-affected area.

Raster scanning
Next, we consider a raster trajectory, adapted to fast writing conditions, typically in a setup configuration where a fast axis is combined with a slow axis. Similarly to the previous experiment, we inscribe a 40 µm × 40 µm square, here using unidirectional parallel lines, separated by a pitch of 0.5 µm, written at a speed of 100 µm s −1 . Given the 120 samples acquired for one line, the average spatial sampling period is of about 0.4 µm (considering the acceleration and deceleration time) and corresponds to one third of the typical modification size. The experimental results are shown in figure 7. Unlike the previous case, a raster approach effectively smooths directional effects originating from the scanning pattern as the beam travels in a single direction. Inscription quality is therefore expected to be higher. Such a method would be preferred for instance for direct waveplate inscription. It is however quite slower, as the sample has to travel back after every move without inscribing.

Towards sub-diffraction resolution patterns
Femtosecond laser exposure being a nonlinear process gives the opportunity to beat the diffraction limit of the focusing optics and to reach feature sizes effectively smaller than the spot size itself. This is achievable when working at pulse energies close to the absorption threshold. However, our in-line measuring setup is still bound to the diffraction limit and cannot resolve features smaller than this limit.
To assess the actual resolution of our inscription technique that unfortunately could not be directly resolved with the DHM, we chose to inscribe a periodic phase-plate down to the limit for which the model still shows satisfying contrast-in our case, 625 or 714 nm (64 and 56 lines in a 40 µm square, respectively). As a reference, the inscription system has a theoretical diffraction limit of 1.3 µm (λ = 1030 nm with an objective lens of NA 0.40). To measure the actual spatial periodicity of the patterns, we implemented a moiré image principle. Technically, we wrote two patterns at a close distance along the optical axis and at a distance smaller than the imaging depth of field. Both patterns are such that they are oriented one relatively to the other with a small and defined angle difference. Superimposed periodic patterns produce moiré fringes that depends on their relative angles. Depending on the angle chosen, these fringes can be spaced at a distance effectively larger than the resolving power of the microscope and are therefore observable. As the angle α is precisely known-it is programmed by numerically tilting the image, the spatial period p of the patterns, can be derived from the period of the moiré fringes Λ, according to Λ≃p/α (for small angle approximation). For both periods, we inscribe two separate sets of patterns with a tilt of either 8 • or 12 • . The same beam path is used in all cases and is only translated to perform the different experiments. The fringes are written orthogonal to the beam direction with a raster pattern. An image of the fringes is fed into the model to generate the modulation signal. The results are shown in figure 8.
The phase maps of the experiments show moiré fringes for different periods. Knowing the angle and the phase map, we can calculate back the actual period of the inscribed structures. In the case of 56 lines, we observe fringes of 12.6, and 8.4 µm for tilting angles of respectively 8 • and 12 • . These values correspond to the exact same value of 1.75 µm, for an expected value of 1.71 µm. In the case of 64 lines, we observe fringes of 7.7, and 6.3 µm for the same angles. These values equate to periods of 1.08, and 1.32 µm, respectively, for an expected value of 1.25 µm.
These experimental observations confirm the ability to inscribe sub-diffraction limit patterns (possibly down to half the measured periods, so here to about 540 nm). However, given the variations between the measurements and the expected values, reaching down to such small patterns is certainly subject to a significant reduction in contrast (i.e. overlap of fringes).

Conclusion and outlook
Femtosecond laser direct-write processes have gained considerably over the last decades, finding applications in numerous fields. However, they remain complex processes strongly dependent on multiple factors, such as beam anisotropy or pulse-front tilt that are intrinsic to any laser writing system, notwithstanding inherent technical difficulties related to effective scanning methods. As a consequence, the reliable and accurate implementation of these processes remain challenging, in particular when it comes to achieving sub-wavelength resolutions or for implementing optimal writing strategies to achieve higher writing throughput. Here, we have shown that (a) using a digital twin to mimic the behaviour of the real system, (b) having sufficient user feedback information (owing to in situ phase-contrast microscopy), and (c) having the ability to perform energy modulation, allows for accelerating the search for optimal exposure pattern and for decreasing static disturbances as well as performing inscription of arbitrary phase maps for an arbitrary laser beam scanning trajectory. Such capability is a necessary step for achieving true sub-wavelength resolution writing of femtosecond lasers index changes in arbitrarily complex patterns. Furthermore, the ability to achieve arbitrary scanning trajectory independent from the final patterns is an essential step towards high yield direct-write processes, necessary for successfully implementing femtosecond laser direct-write processes in realistic industrial settings.