## Abstract

We develop a recurrent neural network framework to model the non-Markovian dynamics exhibited by two-level atoms interacting with the radiation reservoir of a photonic crystal. Despite the strong non-Markovianity of the atomic dynamics induced by the rapid spectral variation in photonic density of states of the photonic reservoir, our recurrent neural network approach is able to capture precise details in the atomic evolution, including the fractional steady-state atomic population inversion and spectral splitting of the atomic transition. We demonstrate the robustness of the recurrent neural network setup against reduced data sets and its effectiveness to deal with systems of increased complexity.

Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

## 1. Introduction

Exploring the dynamics of quantum systems of ever-increasing complexity plays a major role in the development of artificial structures capable of processing information on shorter and shorter time and length scales. Traditionally, the study of quantum systems employs idealised models which consist of small numbers of subsystems and exist within isolation from any external influences. However, such models are not entirely physical as in many if not all situations the interaction with the environment has a major impact on the evolution of the system and cannot be neglected. The field of open quantum systems [1] explores the dynamics of quantum systems coupled to external environments, typically modelled as an ensemble of harmonic oscillators. Using the framework of open quantum systems, we can explore the impact of the environmental degrees of freedom on the dynamics of the system of interest in a formal manner within the master equation approach. This strategy produces a more accurate representation of physical system dynamics, but with the added sophistication comes increased complexity. Even for low-dimensional Hilbert spaces, one is faced with considerable difficulty in unravelling the dynamics induced by the environment. A common approach in reducing this complexity is to deploy the Markovian approximation, which is valid only when the environment relaxes on time scales much shorter than the characteristic time scale of the natural system dynamics. However, in the case of micro-structured photonic systems, such as photonic crystals, fast variations with the frequency of the photonic density of states around the band edges of the photonic band gaps make the Markovian approximation invalid [2,3]. For a two-level atom with excitation energy near to the band edge, the strong coupling between the atomic and the reservoir degrees of freedom results in a dressing of the atomic transitions by photon interactions in which the atom-reservoir attributes become intertwined and results in temporal oscillations and fractional decay of atomic population, spectral splitting and sub-natural line-widths of the atomic transitions [4].

Previous studies utilising machine learning and neural network approaches to quantum system dynamics include studies on implementing machine learning through quantum computing [5–7], applications to quantum many body problems [8,9], phase transitions of condensed matter systems [10,11] and state tomography [12,13]. These past studies clearly demonstrate the promise of applying machine learning techniques within the quantum evolution context.

In our work we extend the formalism developed in Ref. [12] to quantum systems of significantly increased complexity and stronger non-Markovian dynamics. In contrast to previous works that have considered only minor non-Markovian effects introduced by a smooth cut-off in the spectral density, here we deal with non-Markovian effects induced by the square root divergence in the density of states at the boundary of the photonic band gap, which, in practice, can be considered to generate the highest degree of non-Markovianity in the atomic evolution. The structure of Recurrent Neural Networks (RNNs) enables the construction of an RNN master equation [12] which provides the basis of a coherent framework of describing quantum non-Markovian dynamics. Such a master equation has the form of the Markovian Lindblad equation [1], in which the non-Markovian features are introduced through Lindblad operators with a built-in time dependence. This is due to the ability of RNNs to store locally information about previous time evolution steps which can be then used to evaluate future states. In this work, we have developed RNNs that can reproduce the entire time evolution of the two-level atom reduced density matrix for varying values of the detuning of the atomic frequency with respect to the photonic band edge frequency. This outcome demonstrates that the RNN is able to effectively learn the master equation dynamics from a relatively reduced set of training data. The approach introduced here has great potential for application to quantum problems which do not admit analytical solutions and it can be employed to dramatically reduce the time required for high performance computational approaches and to enable exploration of regions of the parameter space currently inaccessible.

This paper is organised as follows: In Sec. 2 we introduce the setup for the RNN system and lay out the various parameters used in its implementation. In Sec. 3, we apply the RNN approach to the single-atom dynamics near the band edge of a photonic crystal and perform quantum state tomography for the isotropic and anisotropic photonic crystal models for fixed and varying values of the atomic frequency detuning with respect to the band edge frequency. Finally, in Sec. 4 we consider the dynamics of the two-atom system with the transition frequencies located in the immediate vicinity of an isotropic photonic band gap model, initially with fixed initial conditions and subsequently with varying initial conditions within the single excitation regime.

## 2. RNN setup

The RNN architecture employed in this study consists of multiple GRU (Gated Recurrent Unit) [14] layers with a final dense layer for the output. The RNN was developed using the Tensorflow 2.3.1 library [15]. Mathematically, we can interpret the internals of the GRU as outputting the time-dependant state $s_t$ through

## 3. RNN modelling of the dynamics of single-atom spontaneous emission in a photonic band gap material

The first system we consider is that of a two-level system with its atomic transition frequency close to a photonic crystal band edge. For the photonic reservoir associated with the photonic band gap material, we use two models for the photonic dispersion, the isotropic and the anisotropic approximations models [2]. The isotropic model of the photonic crystal model assumes that the Brillouin zone boundary is at constant values in $k$ space with $|k| = k_0$, which in practice overestimates the size of the phase space occupied by the Brillouin zone boundary [17]. However, for one dimensional photonic materials this model shows promise and is very tangible within the context of photonic wave guides generated from two-dimensional photonic crystals embedded within three-dimensional structures [18–20]. For convenience, we refer to the relative distance between the transition energy of the two-level system, $\omega _0$, and the top of the band edge energy, $\omega _c$, as the atomic frequency detuning, $\Delta _c = \omega _0 - \omega _c$.

The model we employ to generate the training data is based on the spin-boson model with a Hamiltonian in the rotating-wave approximation given by

In the single excitation sector of the combined Hilbert space we need only consider the following states [21]

where $|{0,1}\rangle_A$ are the ground and excited states for the two level system and $|{0,\lambda}\rangle_B$ are the ground and $\lambda$ mode excited states for the photonic reservoir. As the Hamiltonian conserves excitation number, for any initial state of the form the time evolution is given by and full expressions for the $c_i(t)$ coefficients have been derived in Refs. [2,4].In order to generate the density matrix for the system we simply take the outer product $\rho = |{\phi }\rangle \langle{\phi }|$, and then trace over the reservoir degrees of freedom we obtain the reduced density matrix for the system,

From here on, we use $\tau = \beta t$ as a scaled canonical time of the system and the dimensionless scaled atomic detuning $\Delta = \Delta _c/\beta$, with $\beta$ the effective time scale of the system determined by the photonic crystal configuration [2,4]. Next, using the analytical solutions for the excited state amplitude $c_1(\tau )$ for the (an)isotropic photonic crystal models [2,4], we have generated 3000 training sets for 101 time steps in the scaled time variable $\tau$ for the reduced density matrix, and a further 1000 test samples with varying values for the initial conditions of the system ($c_1(0)$ and $c_0$). Employing this data set, we created a RNN that has 2 GRU layers (with 80D output dimension) and a final dense layer. The number of epochs was set at 400, meaning that each set of training data was used 400 times in refining the RNN’s internal parameter. This process starts with the initial condition $\rho (0)$ and attempts to determine the later time evolution through refinement of the internal parameters by means of gradient descent to minimise the cost function:

where $N$ is the number of samples in the set and $|\cdot |_F$ is the Frobenius norm of the matrix. Here, $\widetilde {\rho }_i$ is the i$^\textrm {th}$ RNN predicted value and $\rho _i$ is the i$^\textrm {th}$ true value from the analytic solution. The minimisation of the $J(\tau )$ functional results directly in the minimisation of the difference between the true values and the RNN’s outputs. Once the RNN is trained on this training data set, we are employing the trained RNN to make predictions made on unseen test data, and subsequently averaged out and normalised the predictions to find the kernel $|\tilde {c}_1(\tau )|^2 = |c_1(\tau )/c_1(0)|^2$. This procedure was performed for fixed detuning $\Delta$ values initially, as such each line in Fig. 2 represents a different instance of the RNN. Clearly, the RNN performs very well and there is very strong agreement between the RNN produced values and the true values for the kernel $|\tilde {c}_1(\tau )|^2$ for both isotropic and anisotropic models of the photonic crystal dispersion (the dotted lines representing the RNN predictions are almost indistinguishable from the thick lines obtained from the analytical solution).In the next step of our investigation, we have allowed the value of the atomic frequency detuning to vary on each sample within the training and testing data sets. As we are introducing a new feature in the input space (the reduced atomic frequency detuning, $\Delta$), we have added an extra GRU layer to the RNN with the same output dimensions as the previous setup. In Fig. 3(a),(c) we have sampled random data sets from the RNNs predictions for unseen detuning values and compared them against the true values for the isotropic and anisotropic cases, respectively. In order to achieve a functional RNN we have appended the detuning value to the inputs of the network, a necessary step as the value of $c_1(\tau )$ is not uniquely determined by its initial condition for varying values of the detuning, $\Delta$. In Fig. 3(b) and (d) we present the variation of the associated cost functions as a function of time. Since the density matrices are obtained through a trace process, the Frobenius norm, used in determining the cost function is bounded

and we note that the average errors in both RNN setups are of the order $10^{-4}$. Despite the clear increase in complexity of the system and in spite of the fact that we have not changed the size of the training set, the RNN predictions are robust enough to accurately reproduce the dynamics of the system. Moreover, the average values of the cost functions for both the anisotropic and isotropic models are of very similar magnitude, despite the additional complexity associated with the isotropic model, whose dynamics is governed by the time evolution of three atom-field dressed states compared to the anisotropic case’s for which only two atom-photon dressed states are present [4]. Our results show that although the isotropic case is more complex, the training data necessarily carries more information and allows for better refinement of the parameter space.We have explored next the influence of the training data set size on the RNN’s predictions for the steady state of the atomic dynamics. In this experiment, we have trained the RNN with an input of the whole data set of the density matrix $\rho (\tau )$ ($\tau \in [0,10]$) and the first quarter ($\tau \in [0,2.5]$) with the added parameter of the atomic frequency detuning value. These were tested over 1000 samples, we note that that there are relatively few values that stray away from the desired steady state solution even as we reduced the amount of data in the system (see Fig. 4). In practice, an additional filtering procedure can be carried out to remove the data points that deviate too far from nearby neighbours and obtain an accurate recreation of the late time dynamics despite a major reduction in the size of the training data set.

## 4. RNN modelling of the dynamics for coupled two-atom systems in a photonic band gap material

In this section, we expand our approach to deal with even stronger non-Markovian dynamics and consider a system consisting of two atoms each interacting with a photonic reservoir associated with a photonic band gap material.

In order to generate our model for the two atom system, we return to Eq. (9) for the single atom. By taking the time derivative, the atomic system dynamics is then described by the following master equation

If we consider two copies of our system together, the joint density operator $\rho _{12} = \rho _1 \otimes \rho _2$ is given by

By adding coupling terms to the Hamiltonian, we can write the two coupled atoms master equation in the form:

The reduced density matrix for the first atom can then be extracted by tracing over the second atom degrees of freedom:

Using the same recurrent neural network design employed for the single-atom case with the detuning value for both atomic frequencies appended to the initial condition and the master equation formalism introduced for the two-atom system, we have generated RNNs to predict the time evolution of the two-atom system. Firstly, we generated data sets with varying values for the detuning in both atomic systems with a fixed initial condition of atom 1 being excited.

It is important to note that the non-Markovian character of the atomic dynamics is strongly amplified in the two-atom system as the single atom dynamics is not only influenced by the coupling of the atomic system to the rapidly varying photonic reservoir, but also by the internal coupling to the secondary atom in which excitations are coherently transferred between the two atoms. The auxiliary atom also exhibits strongly non-Markovian dynamics due to the interaction with the photonic reservoir. In Fig. 5, we present the most and least effective setup predictions for modelling the dynamics from a random selection of testing data. Our results demonstrate a very good agreement between the numerical solutions and the RNNs predictions, proving that even a relatively low-complexity neural network can recapture strongly non-Markovian dynamics.

We have also explored the dependence of the accuracy of the RNN predictions and the size of the RNN and the training data set to find how this would impact the accuracy of our model. The results of this investigation are shown in Fig. 6 and comparing Fig. 6(a) and (c) we note that increasing the size of the GRU layers output dimension without changing the training size actually produced a loss in accuracy. Since we have greatly increased the parameter space, an accurate exploration of its landscape requires a sizable data set. This aspect is further explored in Fig. 6(d) wherein the additional increase in the size of the data training set results in a further improvement in the accuracy.

Also, we note that the Figs. 6(a)-(f) display similar temporal variation as the accuracy decreases at intermediate times before leveling off. We argue that this is due to the fact that in this region the system dynamics transits from being dominated by the photonic reservoir to being dominated by the coherent interactions with the second atom after the photonic reservoir has equilibrated.

Finally, we have expanded the RNN approach to density matrices with arbitrary initial conditions within the single excitation regime. This amounted to increasing the parameter space from the 3D real valued parameter space for $\Delta _{1,2}$ and $\tau$ to the 9D real valued parameter space (three complex amplitudes plus three from $\Delta _{1,2}$ and $\tau$) and, as shown in Fig. 7, this results in less accurate RNN predictions (errors here are of $10^{-2}$ order). The accuracy of the predictions can be easily improved by considering larger training data sets (here we use only 2500 training samples, the same as in the previous example).

## 5. Conclusions

In this paper we have shown that RNNs are a robust tool for modelling the dynamics of highly non-Markovian quantum systems. We demonstrated this for the case of single atom system coupled to the radiation reservoir of a photonic band gap material. The strong variations with the frequency of the photonic reservoir, in particular the divergences associated with the band edge behaviour, lead to highly non-Markovian dynamics which is fully captured within the RNN framework. Furthermore, in the case of a two-atom system coupled to the same strongly non-Markovian photonic reservoir, we show that the RNN approach maintains its efficacy in predicting the complex quantum dynamics. Our work opens the way of exploring the dynamics of systems evolving in even higher dimensional Hilbert spaces, such as $N$-atom models coupled to the radiation reservoir of a photonic crystal for which even strong nonlinear behaviour is predicted [26–28], but for which the parameter space rapidly becomes untenable for conventional numerical approaches.

The RNN framework introduced in this paper is both robust and highly efficient. On a Dell PowerEdge R630 Server (28 cores, 128 GB RAM), the training time required for each setup was about 5 minutes, and the required time to predict over 1000 samples a matter of seconds. Compared to the approximately 30 hours required to generate the two-atom data set, it is orders of magnitude faster. As such it could provide a powerful tool for increasing the speed of such simulations and offer solutions to optimisation problems wherein scanning the whole parameter space by conventional means is not possible. Moreover, it may be possible to deploy the RNN as part of a rapid scanning tool of the parameter space for desired properties that are not accessible either due to cumbersome analysis or emergence through complexity, such as in the case of quantum biological systems where structured vibrational reservoirs can lead to extended coherences in bio-molecular complexes [29]. To this end such tools could prove invaluable in the design and control of quantum systems or in quantifying the relative complexity of different quantum systems.

## Funding

Leverhulme Trust (DS-2017-079); Engineering and Physical Sciences Research Council (EP/L02263X/1, EP/M008576/1, EP/M027791/1).

## Acknowledgments

This work was supported by the Leverhulme Quantum Biology Doctoral Training Centre at the University of Surrey were funded by a Leverhulme Trust training centre grant number DS-2017-079, and the EPSRC (United Kingdom) Strategic Equipment Grant No. EP/L02263X/1 (EP/M008576/1) and EPSRC (United Kingdom) Grant EP/M027791/1 awards to M.F. We acknowledge helpful discussions with the members of the Leverhulme Quantum Biology Doctoral Training Centre.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

The data underlying the findings of this study are available without restriction [30].

## References

**1. **H. P. Breuer and F. Petruccione, * The Theory of Open Quantum Systems* (Oxford University Press, 2002).

**2. **M. Florescu and S. John, “Single-atom switching in photonic crystals,” Phys. Rev. A **64**(3), 033801 (2001). [CrossRef]

**3. **S. John and M. Florescu, “Photonic bandgap materials: towards an all-optical micro-transistor,” J. Opt. A: Pure Appl. Opt. **3**(6), S103–S120 (2001). [CrossRef]

**4. **S. John and T. Quang, “Spontaneous emission near the edge of a photonic band gap,” Phys. Rev. A **50**(2), 1764–1769 (1994). [CrossRef]

**5. **J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, “Quantum machine learning,” Nature **549**(7671), 195–202 (2017). [CrossRef]

**6. **C. Ciliberto, M. Herbster, A. D. Ialongo, M. Pontil, A. Rocchetto, S. Severini, and L. Wossnig, “Quantum machine learning: a classical perspective,” Proc. R. Soc. A **474**(2209), 20170551 (2018). [CrossRef]

**7. **V. Dunjko and H. J. Briegel, “Machine learning and artificial intelligence in the quantum domain: a review of recent progress,” Rep. Prog. Phys. **81**(7), 074001 (2018). [CrossRef]

**8. **G. Carleo and M. Troyer, “Solving the quantum many-body problem with artificial neural networks,” Science **355**(6325), 602–606 (2017). [CrossRef]

**9. **J. Gray, L. Banchi, A. Bayat, and S. Bose, “Machine-learning-assisted many-body entanglement measurement,” Phys. Rev. Lett. **121**(15), 150503 (2018). [CrossRef]

**10. **E. van Nieuwenburg, Y.-H. Liu, and S. Huber, “Learning phase transitions by confusion,” Nat. Phys. **13**(5), 435–439 (2017). [CrossRef]

**11. **J. Carrasquilla and R. G. Melko, “Machine learning phases of matter,” Nat. Phys. **13**(5), 431–434 (2017). [CrossRef]

**12. **L. Banchi, E. Grant, A. Rocchetto, and S. Severini, “Modelling non-Markovian quantum processes with recurrent neural networks,” New J. Phys. **20**(12), 123030 (2018). [CrossRef]

**13. **G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, “Neural-network quantum state tomography,” Nat. Phys. **14**(5), 447–450 (2018). [CrossRef]

**14. **K. Cho, B. van Merrienboer, Ç. Gülçehre, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using RNN encoder-decoder for statistical machine translation,” *CoRR***abs/1406.1078** (2014).

**15. **M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: large-scale machine learning on heterogeneous systems,” (2015). Software available from tensorflow.org.

**16. **D. P. Kingma and J. Ba, “Adam: a method for stochastic optimization,” in 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, Y. Bengio and Y. LeCun, eds. (2015).

**17. **S. John and J. Wang, “Quantum electrodynamics near a photonic band gap: photon bound states and dressed atoms,” Phys. Rev. Lett. **64**(20), 2418–2421 (1990). [CrossRef]

**18. **M. Florescu and S. John, “Resonance fluorescence in photonic band gap waveguide architectures: engineering the vacuum for all-optical switching,” Phys. Rev. A **69**(5), 053810 (2004). [CrossRef]

**19. **M. Florescu, S. Scheel, H. Häffner, H. Lee, D. Strekalov, P. L. Knight, and J. P. Dowling, “Single photons on demand from 3D photonic band-gap structures,” EPL **69**(6), 945–951 (2005). [CrossRef]

**20. **S. Scheel, M. Florescu, H. Häffner, H. Lee, D. V. Strekalov, P. L. Knight, and J. P. Dowling, “Single photons on demand from tunable 3D photonic band-gap structures,” J. Mod. Opt. **54**(2-3), 409–416 (2007). [CrossRef]

**21. **B. M. Garraway, “Nonperturbative decay of an atomic system in a cavity,” Phys. Rev. A **55**(3), 2290–2303 (1997). [CrossRef]

**22. **G. Lindblad, “On the generators of quantum dynamical semigroups,” Commun. Math. Phys. **48**(2), 119–130 (1976). [CrossRef]

**23. **G. McCauley, B. Cruikshank, D. I. Bondar, and K. Jacobs, “Accurate Lindblad-form master equation for weakly damped quantum systems across all regimes,” npj Quantum Inf. **6**(1), 74 (2020). [CrossRef]

**24. **H. Bethe, “Zur Theorie der Metalle,” Eur. Phys. J. A **71**(3-4), 205–226 (1931). [CrossRef]

**25. **M. W. JohnsonM. AminS. Gildert, “Quantum annealing with manufactured spins,” Nature **473**(7346), 194–198 (2011). [CrossRef]

**26. **S. John and T. Quang, “Localization of superradiance near a photonic band gap,” Phys. Rev. Lett. **74**(17), 3419–3422 (1995). [CrossRef]

**27. **S. John and T. Quang, “Photon-hopping conduction and collectively induced transparency in a photonic band gap,” Phys. Rev. A **52**(5), 4083–4088 (1995). [CrossRef]

**28. **S. John and T. Quang, “Quantum optical spin-glass state of impurity two-level atoms in a photonic band gap,” Phys. Rev. Lett. **76**(8), 1320–1323 (1996). [CrossRef]

**29. **A. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. Huelga, and M. Plenio, “The role of non-equilibrium vibrational structures in electronic coherence and recoherence in pigment-protein complexes,” Nat. Phys. **9**, 113 (2013). [CrossRef]

**30. **A. Burgess and M. Florescu, “Dataset for Non-Markovian dynamics in photonic crystals with recurrent neural networks,” figshare (2021), https://doi.org/10.6084/m9.figshare.1473103.