Predictive Modeling of Light–Matter Interaction in One Dimension: A Dynamic Deep Learning Approach

Abstract

The mathematical modeling and the associated numerical simulation of the light–matter interaction (LMI) process are well-known to be quite complicated, particularly for media where several electronic transitions take place under electromagnetic excitation. As a result, numerical simulations of typical LMI processes usually require a high computational cost due to the involvement of a large number of coupled differential equations modeling electron and photon behavior. In this paper, we model the general LMI process involving an electromagnetic interaction medium and optical (light) excitation in one dimension (1D) via the use of a dynamic deep learning algorithm where the neural network coefficients can precisely adapt themselves based on the past values of the coefficients of adjacent layers even under the availability of very limited data. Due to the high computational cost of LMI simulations, simulation data are usually only available for short durations. Our aim here is to implement an adaptive deep learning-based model of the LMI process in 1D based on available temporal data so that the electromagnetic features of LMI simulations can be quickly decrypted by the evolving network coefficients, facilitating self-learning. This enables accurate prediction and acceleration of LMI simulations that can run for much longer durations via the reduction in the cost of computation through the elimination of the requirement for the simultaneous computation and discretization of a large set of coupled differential equations at each simulation step. Our analyses show that the LMI process can be efficiently decrypted using dynamic deep learning with less than 1% relative error (RE), enabling the extension of LMI simulations using simple artificial neural networks.

1. Introduction

Modeling the behavior of electromagnetic waves in free space and also in complex media is well-known to be an involved procedure that requires the concurrent solution of Maxwell’s equations with the Schrödinger wave equation in order to fully capture the wave behavior under time-varying electron dynamics [1,2,3,4]. Various computational methods, such as the Finite-Difference Time-Domain (FDTD) Method, Finite Element Method, Method of Moments, Discontinuous Galerkin Method, etc., have been proposed for solving the Maxwell equations in simple and complex media [5,6,7]. These methods vary in their success in capturing electromagnetic wave behavior in spatial and temporal domains through media that are simple and complex (bi-anisotropic, in-homogenous, nonlinear, etc.).

Under an interaction medium with slowly varying electron dynamics, solving the Maxwell equation is usually sufficient, in which case the electron dynamics can be modeled using the rate equations. The most commonly used method for dealing with complex media in microwave theory, optics, and photonics is the FDTD method. The FDTD method, besides being theoretically simple, is straightforward to implement computationally and can be applied to model wave dynamics in almost all kinds of complex media [5]. However, the FDTD method is limited in its ability to model geometrically complex media or media with rough surfaces and significant shape irregularities. In this study, we will focus on using the FDTD method for modeling the light–matter interaction (LMI) process in nonlinear dispersive media, typically encountered in optics/photonics.

To model the behavior of light, we use the wave equation derived from Maxwell’s equations by substituting the equation for Ampere’s law into the equation for Faraday’s law [5,7]. The dynamics of the electrons are accounted for via the polarization term (also involved in the wave equation) [8,9,10,11,12], which is modeled using the Lorentz dispersion equations [5,6]. The wave equation and the Lorentz dispersion equations are coupled to each other through the polarization term [7]. Under a particle-based treatment of the electrons, the electron dynamics are modeled using the rate equations, helping to identify the population of electrons at each energy level in the medium [5]. The set of all coupled differential equations is modeled based on the FDTD discretization. Importantly, the stability of the numerical solutions relies heavily on the preciseness of the spatial and temporal FDTD discretization. If the discretization step size is too large in the temporal domain, then the solutions can become unstable, meaning that the model can produce non-physical values that jump to infinity [5,6,7]. Hence, for the stability of the numerical solutions, it is very important to choose a sufficiently small temporal step size. Note that the temporal step size must be selected in accordance with the spatial step size. This is because one would be interested in capturing every detail of the propagating light wave through the spatial domain of the involved medium, and for that, the spatial step size should be selected as small (how small depends on the medium properties and the desired resolution of the light wave). Once an adequately small spatial step size is determined, the corresponding temporal step size must satisfy the Courant stability criterion, which states that the temporal step size must be smaller than the spatial step size divided by the speed of light c in free space $∆t<∆x/c$ [5,6,7]. Given this condition, one is left with the conclusion that in order to simulate and predict the light behavior in complex media with high spatial resolution over a long duration, one must select the temporal and spatial step sizes to be very small, which leads to an enormous computational cost [10,11,12].

Here a major problem occurs because the wavelength of optical waves (mid-infrared–visible–ultraviolet–X-ray) is quite small, often less than 1 µm [13,14]. However, many optical phenomena, such as four-wave mixing and parametric amplification, produce a noticeable effect on the macroscale, which is usually much greater than 1 mm. In addition, some optical phenomena, such as interference and laser build-up in a cavity (Figure 1), are only noticeable after relatively long durations (microseconds or milliseconds) [15,16]. In such cases, modeling wave dynamics is a substantial challenge [17,18] which demands the use of supercomputers even under the use of a simple computational algorithm such as FDTD. Therefore, this study aims to come up with a cost-effective alternative solution that extrapolates the solution after a small amount of simulation time.

To achieve this, we resort to deep learning, which can uncover all features of the underlying problem based on a given input-output dataset [19,20,21,22]. Here, we initially define and numerically solve a given problem using the FDTD method under an ultra-fine grid size to ensure stability. Then, for a given input dataset, we achieve an output dataset for a limited simulation duration. We use multilayer perceptrons to form a deep neural network whose coefficients evolve throughout the simulation. After identifying all the coefficients that connect the artificial neurons at each time step, we extrapolate the problem for longer durations by predicting the future values of the network coefficients from their temporal behavior and compare the extrapolated solution with the one that would be attained by running the simulation for the same amount of time for error quantification. Our deep learning model is based on a series of linear-time-invariant (LTI) system expansions that are based on the segmentation of the training data (partitioning of the entire signal data), assuming that the system response changes slowly in time with respect to the period/frequency of the analyzed light wave (which is usually the case in practice) over each of these segments. This is quite accurate because the mathematical model of a multilayer perceptron reduces to multiple LTI system expansions for slowly time-varying systems over segmented training data, where a system can also be approximated as linear [5,6]. Consequently, we use the identified LTI system coefficients at each layer and record their values for all times to extrapolate the wave behavior for longer durations and measure the corresponding extrapolation error by comparing the extrapolated results with the ones obtained via simulation.

Deep learning is based on multilayer neural networks that relate the input data to the output data through a series of coefficient values via inner-product operation [23,24,25,26,27,28,29,30,31,32]. Due to the dramatic increase in the processing power of computers, in the last few decades machine learning has been applied extensively in the field of optics and photonics as it has also been in many other fields [27,28,29]. Over the years, a huge number of papers have focused on machine learning algorithms to uncover hidden physical models in electromagnetics [23,24,25,26,27,28,29]. Notably, research papers that apply machine learning algorithms to predict the nonlinear dynamics of optical systems have received strong attention. These papers have successfully modeled optical systems as robust neural networks that can imitate the behavior of the underlying system with ultimate precision [23,24,25,26,27,28,29,30,31,32]. A striking observation in many of these papers is that the concern is to identify the originally unknown system within the duration of the performed experiments with no extrapolation of the system coefficients for predicting future system behavior. This is problematic because the input signal characteristics often change, leading to a response that differs from the trained model. This unpredictability forms the core challenge in forecasting system behavior, highlighting a crucial battle between traditional machine learning and advanced deep learning algorithms. Deep learning, distinguished from machine learning by its use of multiple layers, offers a solution to this challenge. Each additional layer in a deep learning model extracts more nuanced features of the underlying system, enhancing the model’s predictive capabilities. Predicting the future weights of each layer enables more accurate predictions of subsequent layers’ weights and, consequently, the future samples of an output dataset. This approach requires the adaptive identification of network coefficients over time, segment by segment, a process we refer to as dynamic deep learning. In this study, we will apply dynamic deep learning to both accelerate and predict complex light–matter interaction (LMI) scenarios where conventional prediction methods fall short, enabling accurate and cost-effective simulations in the field of optics and photonics.

2. Methods

2.1. FDTD-Based LMI Model Based on Two Energy Levels

We first investigate the case where photons and electrons interact based on two energy levels with the corresponding frequencies ${\upsilon }_1$ and ${\upsilon }_2$. In this case, there are only two transitions, which are $1\to 2$ and $2\to 1$. Correspondingly, there is only one transition frequency, which is ${\omega }_R=2\pi ({\upsilon }_2-{\upsilon }_1)$. Hence, the polarization density goes into resonance only at this frequency. When the interaction material has a single electronic transition frequency, upon excitation of the interaction medium via an optical beam of electric field intensity E, the corresponding wave, polarization, and electron dynamics are modeled as follows [5,6,7].

• Numerical model for wave dynamics (the wave equation):

$${\nabla }^{2}E-{\mu }_0{\epsilon }_{\mathrm{\infty }}\frac{{\partial }^{2}E}{{\partial t}^2}={\mu }_0{\sigma \frac{\partial E}{\partial t}+\mu }_0\frac{d^{2}P}{{dt}^2}$$

(1)

• Numerical model for polarization dynamics (Lorentz dispersion eq. up to the 3rd order):

$$\frac{d^{2}P}{{dt}^2}+{\gamma }_R\frac{dP}{dt}+{{\omega }_R}^{2}P-\frac{{{\omega }_R}^2{P}^2}{2{(N}_1-N_2)ed}+\frac{{{\omega }_R}^2{P}^3}{{{6{(N}_1-N_2)}^2{e}^{2}d}^2}=\frac{{{(N}_1-N_2)e}^{2}E}{m_e}$$

(2)

• Numerical model for electron dynamics:

$$\frac{d{N}_1}{dt}=-\Gamma N_1+\frac{N_2}{{\tau }_{21}}-\frac{1}{(h/2\pi ){\omega }_R}E\frac{dP}{dt}$$

(3)

$$\frac{d{N}_2}{dt}=\Gamma N_1+\frac{1}{(h/2\pi ){\omega }_R}E\frac{dP}{dt}-\frac{N_2}{{\tau }_{21}}$$

(4)

E: Electric field, μ₀: Free space permeability, ε_∞: Background permittivity, t: Time, h: Planck constant, σ: Conductivity, P: Polarization density, γ_R: Polarization damping rate, ω_R: Angular resonance frequency, N₁: Electron density at the first level, N₂: Electron density at the second level, τ₂₁: Transition lifetime (level 2 to level 1), d: Atom diameter, m_e: Electron mass, e: Elementary charge, and Γ: Electron pump rate.

Here, Equation (1) represents the electric field wave equation for isotropic media that can be easily derived from equations representing Ampere’s Law and Faraday’s Law (also known as Maxwell’s First and Second Equations) [5]. Equation (2) is the famous Lorentz dispersion equation, which governs the temporal behavior of the polarization density of bound electrons under electric field excitation. The Lorentz dispersion equation is analogous to the spring model, which describes the position of a mass that is attached to a spring on which a certain amount of force is exerted. In this case, the electrons are under damped oscillation due to the force that is exerted by the atomic nucleus and the electric field. The damping occurs mainly due to collisions with other electrons. Equations (3) and (4) represent the electron populations at level 1 (ground level) and level 2 (excited level). Since the polarization density in Equation (2) is created by the electron population difference between level 1 and level 2, and since the electric field intensity in Equation (1) depends on the polarization density (which is the source term), Equations (1)–(4) are coupled to each other and must be solved simultaneously. The situation is more complicated in an interaction medium with more energy levels (see next section).

2.2. FDTD-Based LMI Model Based on Three Energy Levels

If there are three energy levels involved in the electron/photon dynamics, then the interaction medium can have six different electron transitions $\left(1\leftrightarrow 2, 1\leftrightarrow 3, 2\leftrightarrow 3\right)$ and three electron transition frequencies (${\upsilon }_3>{\upsilon }_2>{\upsilon }_1,{ \omega }_1=2\pi ({\upsilon }_2-{\upsilon }_1)$, ${\omega }_2=2\pi ({\upsilon }_3-{\upsilon }_1)$, and ${\omega }_3=2\pi ({\upsilon }_3-{\upsilon }_2)$). In this case, the associated LMI process is governed by modified versions of Equations (1)–(4), which are given as follows:

• Numerical model for wave dynamics:

$${\nabla }^{2}E-{\mu }_0{\epsilon }_{\mathrm{\infty }}\frac{{\partial }^{2}E}{{\partial t}^2}={\mu }_0{\sigma \frac{\partial E}{\partial t}+\mu }_0\frac{d^2\left(P_1+P_2+P_3\right)}{{dt}^2}$$

(5)

• Numerical model for polarization dynamics:

$$\frac{d^2{P}_1}{{dt}^2}+{\gamma }_1\frac{d{P}_1}{dt}+{{\omega }_1}^2{P}_1-\frac{{{\omega }_1}^2{P_1}^2}{2{(N}_1-N_2)ed}+\frac{{{\omega }_1}^2{P_1}^3}{{{6{(N}_1-N_2)}^2{e}^{2}d}^2}=\frac{{{(N}_1-N_2)e}^{2}E}{m_e}$$

(6)

$$\frac{d^2{P}_2}{{dt}^2}+{\gamma }_2\frac{d{P}_2}{dt}+{{\omega }_2}^2{P}_2-\frac{{{\omega }_2}^2{P_2}^2}{2{(N}_1-N_3)ed}+\frac{{{\omega }_2}^2{P_2}^3}{{{6{(N}_1-N_3)}^2{e}^{2}d}^2}=\frac{{{(N}_1-N_3)e}^{2}E}{m_e}$$

(7)

$$\frac{d^2{P}_3}{{dt}^2}+{\gamma }_3\frac{d{P}_3}{dt}+{{\omega }_3}^2{P}_3-\frac{{{\omega }_3}^2{P_3}^2}{2{(N}_2-N_3)ed}+\frac{{{\omega }_3}^2{P_3}^3}{{{6{(N}_2-N_3)}^2{e}^{2}d}^2}=\frac{{{(N}_2-N_3)e}^{2}E}{m_e}$$

(8)

• Numerical model for electron dynamics:

$$\frac{d{N}_1}{dt}=-\Gamma N_1+\frac{N_2}{{\tau }_{21}}+\frac{N_3}{{\tau }_{31}}-\frac{1}{(h/2\pi ){\omega }_1}E\frac{d{P}_1}{dt}-\frac{1}{(h/2\pi ){\omega }_2}E\frac{d{P}_2}{dt}, (\Gamma ={\Gamma }_1+{\Gamma }_2)$$

(9)

$$\frac{d{N}_2}{dt}={\Gamma }_1{N}_1+\frac{1}{(h/2\pi ){\omega }_1}E\frac{d{P}_1}{dt}-\frac{N_2}{{\tau }_{21}}+\frac{N_3}{{\tau }_{32}}$$

(10)

$$\frac{d{N}_3}{dt}={\Gamma }_2{N}_1+\frac{1}{(h/2\pi ){\omega }_2}E\frac{d{P}_2}{dt}-\frac{N_3}{{\tau }_{31}}-\frac{N_3}{{\tau }_{32}}$$

(11)

Notice that in this case, the electric field wave equation (Equation (5)) is stimulated by three different polarization terms ($P_1, P_2, P_3$) which are governed by Equations (6)–(8), representing the dynamics of bound charge polarization densities as induced by each electronic transition. The electron populations at each level are modeled via Equations (9)–(11), which are coupled to the Lorentz dispersion equations (Equations (6)–(8)) through the source term and the nonlinear polarization terms. Hence, to solve for the electric field intensity, one must solve all seven equations (Equations (5)–(11)).

2.3. Numerical Discretization

The numerical discretization is carried out in 1D space, considering the case of a single transition based on two-level electron dynamics, using the FDTD method. The discretization for the three-level case is similar and straightforward. Here the electric field and the associated polarization density are discretized in space and time, whereas the electron populations are discretized only in time as it is assumed that the interaction medium is homogenous in density over space. The FDTD discretization of Equations (1)–(4) is carried out as follows [5]:

$$\begin{gathered} x:\mathit{Spatial}\mathit{coordinate},t:\mathrm{Time},E\left(x,t\right)=E\left(i∆x,j∆t\right)\to E\left(i,j\right), i=1,2,\dots ,M_i , j=1,2,\dots ,M_j \\ E\left(x+∆x,t\right)-E\left(x,t\right)\to E\left(i+1,j\right)-E\left(i,j\right), E\left(x,t+∆t\right)-E\left(x,t\right)\to E\left(i,j+1\right)-E\left(i,j\right) \\ \frac{\partial E}{\partial x}=\frac{E\left(i+1,j\right)-E\left(i-1,j\right)}{2∆x},\frac{\partial E}{\partial t}=\frac{E\left(i,j+1\right)-E\left(i,j-1\right)}{2∆t} \\ \frac{{\partial }^{2}E}{{\partial x}^2}=\frac{E\left(i+1,j\right)-2E\left(i,j\right)+E\left(i-1,j\right)}{{∆x}^2},\frac{{\partial }^{2}E}{{\partial t}^2}=\frac{E\left(i,j+1\right)-2E\left(i,j\right)+E\left(i,j-1\right)}{{∆t}^2}, \end{gathered}$$

Notice that we are using the central difference approximation for the first derivatives in time and space and applying the standard formulation for the second derivatives.

• Discretization of the wave equation:

$$\begin{array}{l}E\left(i,j+1\right)=E\left(i+1,j\right)\left\{\frac{2{∆t}^2}{{\mu }_0{∆x}^2\left(\sigma ∆t+2{\epsilon }_{\infty }\left(i,j\right)\right)}\right\}+E\left(i,j\right)\left\{\frac{4{\mu }_0{\epsilon }_{\infty }\left(i,j\right){∆x}^2-2{∆t}^2}{{\mu }_0{∆x}^2\left(\sigma ∆t+2{\epsilon }_{\infty }\left(i,j\right)\right)}\right\}\\ +E\left(i-1,j\right)\left\{\frac{2{∆t}^2}{{\mu }_0{∆x}^2\left(\sigma ∆t+2{\epsilon }_{\infty }\left(i,j\right)\right)}\right\}\\ +E\left(i,j-1\right)\left\{\frac{\left(\sigma ∆t-2{\epsilon }_{\infty }\left(i,j\right)\right)}{\left(\sigma ∆t+2{\epsilon }_{\infty }\left(i,j\right)\right)}\right\}\\ -\frac{2}{\left(\sigma ∆t+2{\epsilon }_{\infty }\left(i,j\right)\right)}\left\{P\left(i,j+1\right)-2P\left(i,j\right)+P\left(i,j-1\right)\right\}\end{array}$$

(12)

• Discretization of the Lorentz dispersion equation:

$$\begin{array}{l}P\left(i,j+1\right)=P\left(i,j-1\right)\left\{\frac{\left(\gamma ∆t-2\right)}{\left(\gamma ∆t+2\right)}\right\}+\frac{2\left(N_1\left(j\right)-N_2\left(j\right)\right)e^2{∆t}^{2}E\left(i,j\right)}{\left(\gamma ∆t+2\right)m}\\ +\frac{2}{\left(\gamma ∆t+2\right)}\{2P\left(i,j\right)-{{\omega }_0}^2{∆t}^{2}P\left(i,j\right)+\frac{{{\omega }_0}^2{{∆t}^{2}P}^2\left(i,j\right)}{\left(N_1\left(j\right)-N_2\left(j\right)\right)ed}\\ -\frac{{{\omega }_0}^2{∆t}^2{P}^3\left(i,j\right)}{{\left(N_1\left(j\right)-N_2\left(j\right)\right)}^2{e}^2{d}^2}\} \end{array}$$

(13)

• Discretization of the rate equations:

$$N_1\left(j+1\right)=N_1\left(j-1\right)+2∆t\left(-\Gamma \left(j\right)N_1\left(j\right)+\frac{N_2\left(j\right)}{{\tau }_{21}}\right)-\frac{2}{h{\omega }_R}E\left(j\right)\left(P\left(j+1\right)-P\left(j-1\right)\right)$$

(14)

$$N_2\left(j+1\right)=N_2\left(j-1\right)+2∆t\left(\Gamma \left(j\right)N_1\left(j\right)-\frac{N_2\left(j\right)}{{\tau }_{21}}\right)+\frac{2}{h{\omega }_R}E\left(j\right)\left(P\left(j+1\right)-P\left(j-1\right)\right)$$

(15)

These discretized equations are solved together with the initial and boundary conditions of a given problem, which will be stated in the analyzed cases in Section 3. For the stability of the numerical solutions, the time step $∆t$ must be chosen to be smaller than the spatial step divided by the speed of light in free space ($∆t<∆x/c)$, which is known as the Courant condition [5].

2.4. Multilayer-Perceptron-Based Neural Networks for Modeling 1D Light-Matter Interaction under the Slowly Varying Envelope Approximation

The FDTD simulations generate an output dataset for each input dataset given by the user. To process these input/output datasets for modeling the underlying system, we consider the simplest neural network that can be mathematically expressed through plain matrix multiplication between the network weights and the input data. In matrix notation, the output data $y$ of a single-layer neural network are expressed in terms of the given input data x (we use brackets for representing network weights and parenthesis for input/output data).

$$\left[\begin{array}{c}y\left(1\right)\\ .\\ y\left(n\right)\end{array}\right]=L_1\left\{\left[\begin{array}{ccc}{w}_1\left[1;1\right]& \cdots & w_1\left[1;K_1\right]\\ ⋮& \ddots & ⋮\\ w_1\left[n;1\right]& \cdots & w_1\left[n;K_1\right]\end{array}\right]\left[\begin{array}{c}x\left(1\right)\\ .\\ x\left(K_1\right)\end{array}\right]\right\}$$

(16)

$$y=L_1\left\{w_{1}x\right\}$$

$$y=\left[y\left(1\right),y\left(2\right),\dots ,y\left(n\right)\right], x=\left[x\left(1\right),x\left(2\right),\dots ,x\left(n\right)\right]$$

where $w_1$ is the coefficient matrix containing the weights that connect/relate the network nodes, and the activation operator $L_1$ can be either a linear or a nonlinear operator that transforms the elements of the result vector of the matrix multiplication to fit within a certain range in accordance with the analyzed problem and the underlying system.

For a given output sample $y\left(n\right)$, Equation (16) can also be expressed via the following summation term:

$$y\left(n\right)=L_1\left\{\sum _{k=1}^{K_1}{w}_1\left[n;k\right]x\left(k\right)\right\}$$

(17)

In system theory, Equation (17) represents the response of a nonlinear time-varying (NTV) system [33,34]. For slowly (time) varying systems, an NTV system can be approximated as a linear time-invariant (LTI) system within narrow time intervals, in which case Equation (17) can be expressed by the time shift between the input and output samples in the form of a convolution operation [34]:

$$y\left(n\right)=L_1\left\{\sum _{k=1}^{K_1}{w}_1\left[k\right]x(n-k)\right\}=L_1\left\{\sum _{k=1}^{K_1}{w}_1\left[n-k\right]x\left(k\right)\right\}, n>K_1$$

(18)

which can be rewritten as a matrix equation ($L_1$ is now an identity operator):

$$\left[\begin{array}{c}y\left(K_1+1\right)\\ .\\ y\left(n\right)\end{array}\right]=L_1\left\{\left[\begin{array}{ccc}{w}_1\left[K_1\right]& \cdots & w_1\left[1\right]\\ ⋮& \ddots & ⋮\\ w_1\left[n-1\right]& \cdots & w_1\left[n-K_1\right]\end{array}\right]\left[\begin{array}{c}x\left(1\right)\\ .\\ x\left(K_1\right)\end{array}\right]\right\}$$

(19)

Equation (19) has much fewer unknowns to solve for (compared to Equation (16)), as the weight (unknown) matrix can now be vectorized as a time shift vector. For this reason, the LTI system approximation provides a great reduction in the computational cost for slowly time-varying systems that are trained gradually over small data segments [35,36,37,38] (see Section 2.5).

In this study, we extend the formulation in Equation (19) for a three-layer neural network, which is considered as a deep learning network with two hidden layers. The use of hidden layers provides more accuracy in determining all features of a given system. In this case, the input–output relationship is modified as follows:

$$\left[\begin{array}{c}y\left(1\right)\\ .\\ y\left(n\right)\end{array}\right]=L_3\left\{\left[\begin{array}{ccc}{w}_3\left[1;1\right]& \cdots & w_3\left[1;K_3\right]\\ ⋮& \ddots & ⋮\\ w_3\left[n;1\right]& \cdots & w_3\left[n;K_3\right]\end{array}\right]L_2\left\{\left[\begin{array}{ccc}{w}_2\left[1;1\right]& \cdots & w_2\left[1;K_2\right]\\ ⋮& \ddots & ⋮\\ w_2\left[K_3;1\right]& \cdots & w_2\left[K_3;K_2\right]\end{array}\right]L_1\left\{\left[\begin{array}{ccc}{w}_1\left[1;1\right]& \cdots & w_1\left[1;K_1\right]\\ ⋮& \ddots & ⋮\\ w_1\left[K_2;1\right]& \cdots & w_1\left[K_2;K_1\right]\end{array}\right]\left[\begin{array}{c}x\left(1\right)\\ .\\ x\left(K_1\right)\end{array}\right]\right\}\right\}\right\}$$

(20)

Based on Equation (20), one can deduce that each output sample is related to the input samples through the following mathematical relation:

$$y\left(n\right)=L_3\left\{\sum _{k=1}^{K_3}{w}_3\left[n;k\right]L_2\left\{\sum _{l=1}^{K_2}{w}_2\left[k;l\right]L_1\left\{\sum _{q=1}^{K_1}{w}_1\left[l;m\right]x\left(m\right)\right\}\right\}\right\}$$

(21)

Just as in the case for a single-layer network, assuming a slowly time-varying system and using small segments of training data, Equation (21) can be expressed in terms of the time shift between the input and output samples over each given training segment (LTI system approximation), which is represented as a serial convolution operation between the network weights and the input data:

$$y\left(n\right)=\sum _{k=1}^{K_3}\sum _{l=1}^{K_2}\sum _{m=1}^{K_1}{w}_3[n-k]{w_2\left[k-l\right]w}_1\left[l-m\right]x\left(m\right)$$

(22)

in which case the activation functions are identity operators. In this study, we will base our approach on Equation (22) to train a 3-layer network for identifying the systems that govern an LMI process.

2.5. Segmentation of a Given Input/Output Dataset for Neural Network Training

In order to evaluate the weights of the three-layer neural network adaptively using many narrow time intervals, one has to divide a given signal into many small segments, each containing M samples that can be used to train the neural network for identifying the local behavior of the underlying system. Each segment should be small compared to the size of the entire signal so that both the linearity and the slowly time-varying system approximation can be made. By computing the coefficients of a 3-layer network at each and every segment and then combining them, one can determine the temporal evolution of the coefficients (thus the network) through the duration of the entire signal, thereby performing dynamic deep learning that facilitates a self-learning network for precise prediction of future output values even under very limited input data. The segmentation can be performed in different ways, as stated below.

• Use a segment window of M samples and a segment period of N (N = M) samples:

Under this segmentation, all segments are adjacent to each other and there is no overlap between them. If the entire signal $S$ is divided into P segments, this is stated as follows:

$$S\left(n;1\right)=\left[s_1,s_2,\dots ,s_M\right], S\left(n;2\right)=\left[s_{M+1},s_{M+2},\dots ,s_{2M}\right],\dots ,S\left(n;P\right)=\left[s_{M\times \left(P-1\right)+1},s_{M\times \left(P-1\right)+2},\dots ,s_{M\times P}\right]$$

(23)

• Use a segment window of M samples and a segment period of N (N $\ne$ M) samples:

For this case, the segments may or may not be overlapping depending on N and M:

$$S\left(n;1\right)=\left[s_1,s_2,\dots ,s_M\right], S\left(n;2\right)=\left[s_{N+1},s_{N+2},\dots ,s_{N+M}\right],\dots ,S\left(n;P\right)=\left[s_{N\times \left(P-1\right)+1},s_{N\times \left(P-1\right)+2},\dots ,s_{N\times \left(P-1\right)+M}\right]$$

(24)

If N < M, the overlap ratio is $(\mathrm{M}-\mathrm{N})/\mathrm{M}$. Whereas if N > M, then the overlap ratio is zero. However, if N is much larger than M, then many samples remain idle, which leads to an inefficient use of data. Hence, if no overlap is desired, one should choose N = M.

• Use a segment window of M samples and a segment period of N (N $\ll$ M) samples:

Here, the overlap ratio is $(\mathrm{M}-\mathrm{N})/\mathrm{M}\approx 1$. This means that each subsequent segment is nearly identical to its neighbor but not completely the same. This is an efficient use of data in cases where the underlying system is strongly time-variant. Whereas if the system is weakly time-variant, this can cause the unnecessary creation of new segments, which would bring little new information about the system (a slow update of the neural network weights would be sufficient) at the cost of drastically increased computation time.

• Use a segment window of M samples and a segment period of N = 1 sample:

In this scenario, each consecutive segment differs only through a single sample shift. This means that subsequent segments totally overlap with each other, except for a single sample such that

$$\begin{gathered} S\left(n;i\right)=\left[s_{M\times \left(i-1\right)},s_{M\times \left(i-1\right)+1},\dots ,s_{M\times i-1}\right],S\left(n;i+1\right)=\left[s_{M\times \left(i-1\right)+1},s_{M\times \left(i-1\right)+2},\dots ,s_{\left(M\times i\right)}\right] \\ i=1,2,\dots P \left(i:Segment number\right) \end{gathered}$$

(25)

This form of segmentation (Figure 2) enables a very efficient use of data under fast time-variation of the system coefficients. Despite our assumption of a slowly time-varying system, we will make use of this segmentation for a sample-by-sample computation of the neural network coefficients in order to see their temporal evolution in full resolution.

Under such sample-by-sample segmentation, to predict the upcoming samples of the output vector y, one has to solve for the set of coefficients at each layer and for each estimated sample. This means that to predict N future samples, the coefficient vectors $w_1,w_2,w_3$ must be solved based on Equation (22), which can be expanded as three consecutive matrix multiplications:

$$\left[\begin{array}{c}y\left(M+1\right)\\ .\\ y\left(M+N\right)\end{array}\right]=\left[\begin{array}{ccc}{w}_3\left[1\right]& \cdots & w_3\left[K_3\right]\\ ⋮& \ddots & ⋮\\ w_3\left[N\right]& \cdots & w_3\left[N-K_3+1\right]\end{array}\right]\left[\begin{array}{c}z\left(M\right)\\ .\\ z\left(M-K_3+1\right)\end{array}\right]$$

(26)

$$\left[\begin{array}{c}z\left(M-K_3+1\right)\\ .\\ z\left(M\right)\end{array}\right]=\left[\begin{array}{ccc}{w}_2\left[1\right]& \cdots & w_2\left[K_2\right]\\ ⋮& \ddots & ⋮\\ w_2\left[K_3\right]& \cdots & w_2\left[K_3-K_2\right]\end{array}\right]\left[\begin{array}{c}q\left(M-K_3\right)\\ .\\ q\left(M-K_3-K_2+1\right)\end{array}\right]$$

(27)

$$\left[\begin{array}{c}q\left(M-K_3-K_2+1\right)\\ .\\ q\left(M-K_3\right)\end{array}\right]=\left[\begin{array}{ccc}{w}_1\left[1\right]& \cdots & w_1\left[K_1\right]\\ ⋮& \ddots & ⋮\\ w_1\left[K_2\right]& \cdots & w_1\left[K_2-K_1\right]\end{array}\right]\left[\begin{array}{c}x\left(M-K_3-K_2\right)\\ .\\ x\left(M-K_3-K_2-K_1+1\right)\end{array}\right]$$

(28)

where x is the input vector and {q, z} are intermediate outputs that ultimately yield the output vector y. Simplified forms of Equations (26)–(29) will be referred to for carrying out the dynamic coefficient (system) identification procedures in the next section based on the least-squares method, using sample-by-sample segmentation under the use of linear operators. Note that Equations (26)–(28) imply $M\ge \left(K_3+K_2+K_1\right), K_3>K_2, \mathrm{a}\mathrm{n}\mathrm{d} K_2>K_1$.

2.6. Training of the Multilayer Perceptron over Each Data Segment Based on Output Data

The system coefficients can be evaluated solely based on past samples of the output data $y$. Hence, one does not necessarily need the input data to perform a prediction. Also, one can simply use Equation (22) over each training segment where nonlinear operators are not required due to the limited size of the involved segment. By re-labeling the coefficients as a, b, and c, Equation (22) can be rewritten as follows (note that x[m] is replaced with y[m] as past samples of the output data are considered input signal samples in the case of prediction):

$$y\left(n\right)=\left\{\sum _{k=n-K_3}^{n-1}a[n-k]\left\{\sum _{l=k-K_2}^{k-1}b\left[k-l\right]\left\{\sum _{m=l-K_1}^{l-1}c\left[l-m\right]y\left(m\right)\right\}\right\}\right\}, n>k$$

(29)

Using the following operational tricks, Equation (29) can be converted into Equation (36):

$$y\left(n\right)=\sum _{k=n-K_3}^{n-1}a[n-k]z\left(k\right)=\sum _{k=n-K_3}^{n-1}a\left[k\right]z(n-k)$$

(30)

$$y\left(n\right)=\left\{\sum _{k=1}^{K_3}a[k]\left\{\sum _{l=n-k-1}^{n-k-K_2}b\left[n-k-l\right]\left\{\sum _{m=1}^{K_1}c\left[m\right]y(l-m)\right\}\right\}\right\}$$

(31)

$$l\leftrightarrow n-k-l$$

$$y\left(n\right)=\left\{\sum _{k=1}^{K_3}a[k]\left\{\sum _{l=1}^{K_2}b\left[l\right]\left\{\sum _{m=1}^{K_1}c\left[m\right]y(n-k-l-m)\right\}\right\}\right\}$$

(32)

$$k\leftrightarrow k-l$$

$$y\left(n\right)=\left\{\sum _{k=l+1}^{l+K_3}a[k-l]\left\{\sum _{l=1}^{K_2}b\left[l\right]\left\{\sum _{m=1}^{K_1}c\left[m\right]y(n-k-m)\right\}\right\}\right\}$$

(33)

$$l\leftrightarrow l-m$$

$$y\left(n\right)=\left\{\sum _{k=l-m+1}^{l-m+K_3}a[k-l+m]\left\{\sum _{l=m+1}^{m+K_2}b\left[l-m\right]\left\{\sum _{m=1}^{K_1}c\left[m\right]y(n-k-m)\right\}\right\}\right\}$$

(34)

$$k\leftrightarrow k-m$$

$$y\left(n\right)=\left\{\sum _{k=l+1}^{l+K_3}a[k-l]\left\{\sum _{l=m+1}^{m+K_2}b\left[l-m\right]\left\{\sum _{m=1}^{K_1}c\left[m\right]y(n-k)\right\}\right\}\right\}$$

(35)

$$y\left(n\right)=\sum _{k=l+1}^{l+K_3}\sum _{l=m+1}^{m+K_2}\sum _{m=1}^{K_1}c\left[m\right]b\left[l-m\right]a\left[k-l\right]y(n-k)$$

(36)

The good thing about Equation (36) is that it can be nicely decoupled into the following three equations:

$$y\left(n\right)=\sum _{k=1}^{K_3}a\left[k\right]y(n-k)$$

(37)

$$a\left[k\right]=\sum _{l=1}^{K_2}b\left[l\right]a\left[k-l\right]$$

(38)

$$b\left[l\right]=\sum _{m=1}^{K_1}c\left[m\right]b\left[l-m\right]$$

(39)

where the output signal can be directly expressed in terms of its past values in a single equation.

Equations (37)–(39) can be rewritten using matrix equations that describe the prediction of N upcoming samples based on the past M samples (segmental training data):

$$\left[\begin{array}{c}y\left(M+1\right)\\ .\\ y\left(M+N\right)\end{array}\right]=\left[\begin{array}{ccc}y\left[M\right]& \cdots & y\left[M+1-K_3\right]\\ ⋮& \ddots & ⋮\\ y\left[M+N-1\right]& \cdots & y\left[M+N-K_3\right]\end{array}\right]\left[\begin{array}{c}a\left[1\right]\\ .\\ a\left[K_3\right]\end{array}\right]$$

(40)

$$\left[\begin{array}{c}a[K_2+1]\\ .\\ a{[K}_3]\end{array}\right]=\left[\begin{array}{ccc}a\left[K_2\right]& \cdots & a\left[1\right]\\ ⋮& \ddots & ⋮\\ a\left[K_3-1\right]& \cdots & a\left[K_3-K_2\right]\end{array}\right]\left[\begin{array}{c}b\left[1\right]\\ .\\ b\left[K_2\right]\end{array}\right]$$

(41)

$$\left[\begin{array}{c}b[K_1+1]\\ .\\ b{[K}_2]\end{array}\right]=\left[\begin{array}{ccc}b\left[K_1\right]& \cdots & b\left[1\right]\\ ⋮& \ddots & ⋮\\ b\left[K_2-1\right]& \cdots & b\left[K_2-K_1\right]\end{array}\right]\left[\begin{array}{c}c\left[1\right]\\ .\\ c\left[K_1\right]\end{array}\right]$$

(42)

Using Equation (40), one can solve for the set of $K_3$ coefficients for the first layer via the least-squares method. The same goes for the second- and third-layer coefficients. The coefficient matrices always depend on the past values of the coefficients of the previous layer ($K_2$ coefficients for the second layer and $K_1$ coefficients for the first layer). Once all the coefficients for each layer are solved, future samples of y(n) can be predicted using sequential matrix updates and multiplication. However, Equations (40)–(42) do not yield the time-dependent values of the coefficients, which prevents the prediction of future values of the coefficients and hinders the adaptive prediction of future output samples. For a self-learning or self-adaptive neural network, the coefficients must evolve with each and every new sample so that the future values of the coefficients can also be predicted from their time series, enabling more accurate output signal prediction under the limited availability of input data. This enables a dynamic evolution of multilayer perceptrons through segmental training where each segment differs only via a single sample shift so that the coefficients are instantly adapted for an “infinitesimal” (unit sample) change. To attain the matrices that contain the time variation (past values) of the coefficients, every set of coefficients needs to be solved per new sample. Therefore, over each segment, the training of the network should be carried out as follows:

$$\left[\begin{array}{c}y\left[n-T-M\right]\\ .\\ y\left[n-T-1\right]\end{array}\right]=\left[\begin{array}{ccc}y\left[n-T-M-1\right]& \cdots & y\left[n-T-M-K_3\right]\\ ⋮& \ddots & ⋮\\ y\left[n-T-2\right]& \cdots & y\left[n-T-K_3-1\right]\end{array}\right]\left[\begin{array}{c}{a}_1\left[n-T-1\right]\\ .\\ a_{K_3}\left[n-T-1\right]\end{array}\right]$$

(43)

$$\left[\begin{array}{c}{a}_1\left[n-T-1\right]\\ .\\ a_{K_3}\left[n-T-1\right]\end{array}\right]=\left[\begin{array}{ccc}{a}_1\left[n-T-2\right]& \cdots & a_1\left[n-T-K_2-1\right]\\ ⋮& \ddots & ⋮\\ a_{K_3}\left[n-T-2\right]& \cdots & a_{K_3}\left[n-T-K_2-1\right]\end{array}\right]\left[\begin{array}{c}{b}_1\left[n-T-1\right]\\ .\\ b_{K_2}\left[n-T-1\right]\end{array}\right]$$

(44)

$$\left[\begin{array}{c}{b}_1\left[n-T-1\right]\\ .\\ b_{K_2}\left[n-T-1\right]\end{array}\right]=\left[\begin{array}{ccc}{b}_1\left[n-T-2\right]& \cdots & b_1\left[n-T-K_1-1\right]\\ ⋮& \ddots & ⋮\\ b_{K_2}\left[n-T-2\right]& \cdots & b_{K_2}\left[n-T-K_1-1\right]\end{array}\right]\left[\begin{array}{c}{c}_1\left[n-T-1\right]\\ .\\ c_{K_1}\left[n-T-1\right]\end{array}\right]$$

(45)

$$n:Sample number, M:Segment size, T:Time shift, \mathrm{T}=1,2,3,4,\dots ,\mathrm{L}-\mathrm{M}$$

Using Equations (43)–(45), one can perform an adaptive prediction at every instantaneous sample of the output signal based on the 3-layer network. This requires the instantaneous solution of each set of coefficients through the least-squares method, as follows:

$$y=Y^{T}a, a=A^{T}b, {b=B}^{T}c, A, B, Y:Coefficient matrices$$

(46)

$$a={{(Y}^{T}Y)}^{-1}{Y}^{T}y, b={{(A}^{T}A)}^{-1}{A}^{T}a, c={{(B}^{T}B)}^{-1}{B}^{T}B$$

(47)

Therefore, once the coefficients over a certain segment are identified, the future samples of the output signal y can be predicted via Equation (48):

$$y=Y^T{A}^T{B}^{T}c={\left(BAY\right)}^{T}c$$

(48)

Note that the coefficient matrices are constructed from past samples of the coefficients themselves. Hence, once a certain coefficient vector (a, b) is solved, its coefficient matrix (A, B) is known, which enables the prediction in Equation (48).

3. Results

3.1. Frequency Generation via Optical Wave Mixing in a Micro-Resonator

3.1.1. Difference Frequency Generation Process

The first optical process that we model using a three-layer network is the difference frequency generation process, which is a nonlinear process that occurs when two high-amplitude waves $E_1$ and $E_2$ are simultaneously fed into a cavity. Here, the cavity is composed of two perfectly reflective walls and a Gallium–Arsenide interaction medium that is inserted between the cavity walls. The excited waves are given as follows:

$$\begin{gathered} E_1\left(x=0 \mathsf{\mu }\mathrm{m},t=0 \mathrm{s}\right)=4\times {10}^8\times \mathit{sin}\left(2\pi \left(2.5\times {10}^{14}\right)t\right) \mathrm{V}/\mathrm{m} (Infrared wave at 250 \mathrm{T}\mathrm{H}\mathrm{z}) \\ {E}_2\left(x=0 \mathsf{\mu }\mathrm{m},t=0 \mathrm{s}\right)=3\times {10}^8\times \mathit{sin}\left(2\pi \left(1.5\times {10}^{14}\right)t\right) \mathrm{V}/\mathrm{m} (Infrared wave at 150 \mathrm{T}\mathrm{H}\mathrm{z}) \end{gathered}$$

Hence, the total excitation wave at t = 0 s is given as follows:

$$\begin{gathered} E_{total}\left(x=0 \mathsf{\mu }\mathrm{m},t=0 \mathrm{s}\right) \\ =4\times {10}^8\times \mathit{sin}\left(2\pi \left(2.5\times {10}^{14}\right)t\right)\mathrm{V}/\mathrm{m}+3\times {10}^8\times \mathit{sin}\left(2\pi \left(1.5\times {10}^{14}\right)t\right) \mathrm{V}/\mathrm{m} \end{gathered}$$

The simulation setting and the properties of the interaction medium are given as follows:

$$\begin{gathered} Spatial and temporal discretization: 0\le x\le 5 \mathsf{\mu }\mathrm{m}, 0\le t\le 5 \mathrm{p}\mathrm{s}, ∆x=10 \mathrm{n}\mathrm{m}, ∆t=0.025 \mathrm{f}\mathrm{s} \\ Interaction medium resonance frequency: f_r=3.45\times {10}^{14} \mathrm{H}\mathrm{z}, Conductivity:\mathsf{\sigma }={10}^{-4} \mathrm{S}/\mathrm{m} \\ Polarization damping rate of the interaction medium: \mathsf{\gamma }=1\times {10}^{10} \mathrm{H}\mathrm{z} \\ Background permittivity: \left({\epsilon }_{\mathrm{\infty }}\right)=13{\epsilon }_0, Atomic diameter: d=0.3 \mathrm{n}\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{s} \\ {\epsilon }_0:Free space permittivity \\ Cavity range:0 \mathsf{\mu }\mathrm{m}<\mathrm{x}<5 \mathsf{\mu }\mathrm{m}, Electron density: \mathrm{N}=1.76\times {10}^{29}/{\mathrm{m}}^3 \end{gathered}$$

As a result of the high optical intensity and the inherent nonlinearity of Gallium–Arsenide atoms, after a certain time, new harmonics will be generated at 100 THz (difference harmonic, mid-infrared) and 400 THz (sum harmonic, red) due to the nonlinear (cross) polarization terms in Equation (2). At the end of the simulation, the final spectrum of the intracavity wave (magnitude of the Fourier transform) is as given in Figure 3.

The modeling of the wave mixing process is computationally costly, and the cost of computation grows exponentially with increasing simulation duration. The efficiency of our three-layer neural network is tested in this case by predicting the temporal behavior of the electric field amplitude of the difference frequency (100 THz) from t = 5 ps to t = 8.5 ps.

The formulations in Equations (43)–(45) are applied with n = 42,500, M = 25,000, $K_3=$ 4, $K_2=$ 4, and $K_1=$ 4 (each layer has four coefficients). Hence, the available training data consist of M = 25,000 samples, and the number of samples to be predicted is 17,500 (n-M). The resulting prediction is shown in Figure 4. In this case, the prediction is highly accurate, with less than 0.08% relative error (RE). A comparison with the theoretical results is also given in Figure 5 regarding the electric field amplitude, based on the given formulation in [39]. The comparison is performed by considering the simulation duration as the total number of roundtrips multiplied by two times the cavity length times the square root of the background permittivity of the medium, divided by the speed of light.

The time variation of the set of four coefficients ($a_1,a_2,a_3,a_4$) for the first layer of the network is given in Figure 6. Notice that the coefficients are nearly constant as the change in percentage around their mean value is quite insignificant.

The coefficients of the second layer are shown in Figure 7. Compared to the coefficients of the first layer, the percentage of change between the minimum and maximum values of these coefficients is much greater. However, the coefficients of the second layer tend to converge to a certain value after 5 picoseconds, which indicates that the coefficients of the first layer will continue to make the same tiny oscillations around their mean values, indicating that a reliable prediction in the long term is feasible.

Finally, the coefficients of the third layer are given in Figure 8. These coefficients are used to predict whether the convergent behavior of the coefficients of the second layer will persist. There is also a tendency for convergence in this layer, which is easier to spot for the second, third, and fourth coefficients, whereas for the first coefficient, convergence is not that critical as the percentage of change is smaller. Therefore, the third-layer coefficients indicate that the convergence behavior of the second-layer coefficients will persist, which further indicates that the extremely small percentage of change for the first-layer coefficients will continue, enabling one to treat them as constant variables over time and perform a prediction for a duration that is much greater than t = 5 ps.

3.1.2. Sum Frequency Generation Process

As a result of the wave mixing process, the time evolution of the generated electric field of the sum harmonic at 400 THz is illustrated in Figure 9 for both the available data (0 < t < 5 ps) and the predicted data (5 ps < t < 8.5 ps). Expectedly, the electric field of the sum harmonic frequency gradually increases from 0 to ${10}^9$V/m in 8.5 ps. The available data suggest that the field amplitude reaches ${4\times 10}^7$V/m in 5 picoseconds, whereas the predicted data show that this amplitude will reach ${10}^9$V/m in the next 3.5 picoseconds. Due to the nonlinearity of the problem, the time it takes to complete the simulation in 8.5 ps is 3.7 times greater than the time it takes for the simulation to be completed in 5 ps. With the aid of the computed neural network coefficients, the entire process is simulated at a duration that is 100 times smaller. The resulting RE of the prediction is below 0.05% (Equations (43)–(45) are used with n = 42,500, M = 25,000, $K_3=$ 4, $K_2=$ 4, and $K_1=$ 4). A comparison with the theoretical results is given in Figure 10 regarding the electric field amplitude, based on the given formulation in [39] using the same trick as discussed in Section 3.1.1.

The time variation of the set of four coefficients ($a_1,a_2,a_3,a_4$) for the first layer of the network is given in Figure 11. Just as in the case of difference frequency generation, the coefficients are nearly constant in time, whereas the coefficients of the second layer (see Figure 12) are not constant in time. However, the coefficients of the second layer converge to a constant value over time, indicating that the coefficients of the first layer will behave similarly in the future.

Although it is not as evident as it is for the second-layer coefficients, the coefficients of the third layer also display a convergent behavior (Figure 13), supporting the stability of the second-layer coefficients around their converged value in the long term. The coefficients of the third layer can also be used to verify the convergent behavior of the second-layer coefficients through insertion in the expansion formula in Equation (39).

3.2. Optical Parametric Amplification

Next, we model the optical parametric amplification (OPA) process using our three-layer network model. Here, as opposed to the two previous cases (frequency generation), only one of the waves has a high amplitude, while the other usually has a low amplitude. Through the wave mixing process, over time, the low-amplitude wave draws energy from the cavity that is energized by the high-amplitude wave and becomes amplified while the high-amplitude wave is attenuated. In this case, we consider the amplification of an infrared signal wave at 150 THz using a high-amplitude pump wave at 250 THz. The cavity configuration remains the same as in the previous section. The signal and the pump waves are given as follows (the “signal” and “pump” designations are specific to OPA):

$$\begin{gathered} E_1\left(x=0 \mathsf{\mu }\mathrm{m},t=0 \mathrm{s}\right)=4\times {10}^8\times \sin \left(2\pi \left(2.5\times {10}^{14}\right)t\right) \mathrm{V}/\mathrm{m} (pump wave) \\ {E}_2\left(x=0 \mathsf{\mu }\mathrm{m},t=0 \mathrm{s}\right)=1\times \sin \left(2\pi \left(1.5\times {10}^{14}\right)t\right) \mathrm{V}/\mathrm{m} (signal wave to be amplified) \end{gathered}$$

Therefore, the total excitation wave at t = 0 s is given as follows:

$$E_{total}\left(x=0 \mathsf{\mu }\mathrm{m},t=0 \mathrm{s}\right)=4\times {10}^8\times \sin \left(2\pi \left(2.5\times {10}^{14}\right)t\right)\mathrm{V}/\mathrm{m}+1\times \sin \left(2\pi \left(1.5\times {10}^{14}\right)t\right) \mathrm{V}/\mathrm{m}$$

The time evolution of the electric field is given in Figure 14 for both the computed and the predicted data. The RE in prediction accuracy was found to be 0.77% in this case. As before, the prediction is performed from t = 5 ps to t = 8.5 ps using the computed data from t = 0 s to t = 5 ps. Once again, the formulations in Equations (43)–(45) are applied with n = 42,500, M = 25,000, ${\mathrm{K}}_3=$ 4, ${\mathrm{K}}_2=$ 4, and ${\mathrm{K}}_1=$ 4. The simulation time using the physical model was 141 s, whereas the computation time is 2.2 s when the derived neural network model is used.

The coefficients of the first layer of the neural network are given in Figure 15. Despite persistent fluctuations, the values of the coefficients are nearly constant over time. The coefficients of the second layer (Figure 16) indicate a convergent behavior, hinting that the temporal behavior of the process is stabilized. The third-layer coefficients (Figure 17) are observed to fluctuate within a narrow range and are stable in their temporal behavior, which hints for a lasting convergence of the second-layer coefficients, implying that the first-layer coefficients will continue to remain (nearly) constant over time.

3.3. Time Evolution of Optical Wave Amplitude in Lossy Micro-Resonators

Optical micro-resonators (or micro-cavities) are often used to confine light waves for sufficiently long durations in order to enable the effective occurrence of certain processes such as parametric amplification, supercontinuum generation, etc. For this reason, the degree of loss in a micro-resonator is of concern as it leads to a decrease in the photon lifetime. In this part, we investigate a lossy optical micro-resonator with an increased interaction medium conductivity and a greater polarization damping rate and predict the overall decay rate for the electric field amplitude within the micro-resonator, for which an analytical expression is not available in the presence of nonlinearity. Here, the total excitation wave at t = 0 s is given as follows:

$$E_{total}\left(x=0 \mathsf{\mu }\mathrm{m},t=0 \mathrm{s}\right)=1\times {10}^6\times \mathit{sin}\left(2\pi \left(2.5\times {10}^{14}\right)t\right)\mathrm{V}/\mathrm{m} (Infrared wave at 250 \mathrm{T}\mathrm{H}\mathrm{z})$$

The simulation setting and the properties of the interaction medium are given as follows:

$$\begin{gathered} Spatial and temporal discretization: 0\le x\le 5 \mathsf{\mu }\mathrm{m}, 0\le t\le 5 \mathrm{p}\mathrm{s}, ∆x=10 \mathrm{n}\mathrm{m}, ∆t=0.025 \mathrm{f}\mathrm{s} \\ Interaction medium resonance frequency: f_r=3.45\times {10}^{14} \mathrm{H}\mathrm{z} \\ Polarization damping rate of the interaction medium: \gamma =1\times {10}^{11} \mathrm{H}\mathrm{z} \\ Interaction medium conductivity: \sigma ={10}^{-2} \mathrm{S}/\mathrm{m} \\ Background permittivity: \left({\epsilon }_{\mathrm{\infty }}\right)=13{\epsilon }_0, Atomic diameter: d=0.3 \mathrm{n}\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{s} \\ Cavity range:0 \mathsf{\mu }\mathrm{m}<x<5 \mathsf{\mu }\mathrm{m}, Electron density: N=1.76\times {10}^{29}/{\mathrm{m}}^3 \end{gathered}$$

For the sake of simplicity, the walls of the micro-cavity are assumed to be perfectly reflective, so the cavity loss is solely due to the dielectric and conduction losses in the interaction medium. The resulting computed and predicted data are shown in Figure 18. Using a three-layer neural network (Equations (43)–(45) are used with n = 65,000, M = 25,000, $K_3=$ 16, $K_2=$ 16, and $K_1=$ 16) with 16 coefficients at each layer, it is successfully (RE: 0.94%) predicted that the electric field amplitude in the micro-cavity will decrease by 12% within 13 picoseconds, which enables an estimation of the total decay rate in the cavity.

The time variation of the first-, second-, and third-layer coefficients are given in Figure 19, Figure 20 and Figure 21, respectively. Compared to the previously analyzed cases, the coefficients for each layer are more fluctuating due to the dominance of the interference effect in the cavity.

4. Conclusions

The light–matter interaction (LMI) process can be accurately decrypted using a three-layer neural network that can be trained during the simulation of a given LMI process that is based on the actual physical model. The network coefficients are dynamically evaluated over many consecutive segments of the available simulation data (each consecutive segment differs only by a sample shift), and their corresponding time evolution can be attained by recording their values for each segment. Based on the time evolution of the network coefficients of the second and third layers, the upcoming values of the first layer can be predicted, which can be used to extrapolate the future values of the output signal. The extrapolation can be made simply through a convolution operation between the network weights and the past values of the output data since the training of the network is performed in a piecewise manner through segmentation, where a linearization approximation can be made. Our computations show that the extrapolation error is quite low for the predicted part of the electric field (light) amplitude, which suggests that the use of a well-trained adaptive neural network can drastically reduce the computation time for the optical electric field and can be used to accelerate long simulations in optics and photonics for modeling scenarios where ultra-fine spatial and/or temporal grid-size is required.