# Highly efficient storage of 25-dimensional photonic qudit in a cold-atom-
based quantum memory

Ming-Xin Dong  Key Laboratory of Quantum Information, University of Science
and Technology of China, Hefei, Anhui 230026, China.  Synergetic Innovation
Center of Quantum Information and Quantum Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China.  School of Physics and
Materials Engineering, Hefei Normal University, Hefei, Anhui 230601, China.
Wei-Hang Zhang  Key Laboratory of Quantum Information, University of Science
and Technology of China, Hefei, Anhui 230026, China.  Synergetic Innovation
Center of Quantum Information and Quantum Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China.  Lei Zeng  Key Laboratory of
Quantum Information, University of Science and Technology of China, Hefei,
Anhui 230026, China.  Synergetic Innovation Center of Quantum Information and
Quantum Physics, University of Science and Technology of China, Hefei, Anhui
230026, China.  Ying-Hao Ye  Key Laboratory of Quantum Information, University
of Science and Technology of China, Hefei, Anhui 230026, China.  Synergetic
Innovation Center of Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei, Anhui 230026, China.  Da-Chuang Li
<EMAIL_ADDRESS>School of Physics and Materials Engineering, Hefei
Normal University, Hefei, Anhui 230601, China.  Guang-Can Guo  Key Laboratory
of Quantum Information, University of Science and Technology of China, Hefei,
Anhui 230026, China.  Synergetic Innovation Center of Quantum Information and
Quantum Physics, University of Science and Technology of China, Hefei, Anhui
230026, China.  Dong-Sheng Ding<EMAIL_ADDRESS>Key Laboratory of Quantum
Information, University of Science and Technology of China, Hefei, Anhui
230026, China.  Synergetic Innovation Center of Quantum Information and
Quantum Physics, University of Science and Technology of China, Hefei, Anhui
230026, China.  Bao-Sen Shi<EMAIL_ADDRESS>Key Laboratory of Quantum
Information, University of Science and Technology of China, Hefei, Anhui
230026, China.  Synergetic Innovation Center of Quantum Information and
Quantum Physics, University of Science and Technology of China, Hefei, Anhui
230026, China.

###### Abstract

Building an efficient quantum memory in high-dimensional Hilbert spaces is one
of the fundamental requirements for establishing high-dimensional quantum
repeaters, where it offers many advantages over two-dimensional quantum
systems, such as a larger information capacity and enhanced noise resilience.
To date, there have been no reports about how to achieve an efficient high-
dimensional quantum memory. Here, we experimentally realize a quantum memory
that is operational in Hilbert spaces of up to 25 dimensions with a storage
efficiency of close to 60%. The proposed approach exploits the spatial-mode-
independent interaction between atoms and photons which are encoded in
transverse size-invariant orbital angular momentum modes. In particular, our
memory features uniform storage efficiency and low cross-talk disturbance for
25 individual spatial modes of photons, thus allowing storing arbitrary qudit
states programmed from 25 eigenstates within the high-dimensional Hilbert
spaces, and eventually contributing to the storage of a 25-dimensional qudit
state. These results would have great prospects for the implementation of
long-distance high-dimensional quantum networks and quantum information
processing.

Introduction. Quantum memories lvovsky2009optical ; sangouard2011quantum  that
enable quantum state storage and its on-demand retrieval are essential
requirements for quantum-repeater-based quantum communication networks
briegel1998quantum ; kimble2008quantum  and scalable quantum computation
kok2007linear . The storage efficiency exceeding the 50% threshold
vernaz2018highly ; wang2019efficient ; guo2019high  is necessary for practical
applications due to the fundamental requirements of beating the quantum no-
cloning limit without post-selection grosshans2001quantum  or realizing error
correction in linear optical quantum computation varnava2006loss . Although
quantum memory has been widely demonstrated in conventional two-dimensional
(or qubit) quantum systems, it is highly desirable to realize a high-
dimensional quantum memory since manipulating a photon in a high-dimensional
Hilbert space, i.e., qudit, provides many advantages over the qubit systems in
terms of practical quantum information processing. For example, qudits enable
networks to carry more information and increase their channel capacity via
superdense coding in quantum communication erhard2018twisted ;
erhard2020advances ; krenn2014generation ; for quantum cryptography, it has
been shown that qudits can provide a more secure flux of information against
eavesdroppers bechmann2000quantum ; bechmann2000quantumprl ; cerf2002security
; walborn2006quantum  since the upper bound of limited cloning fidelity, given
by $F_{{\rm clon}}^{d}=1/2+1/(d+1)$, scales inversely with the dimension
erhard2018twisted , and they also feature a better resilience to noise
sheridan2010security ; ecker2019overcoming . Moreover, qudit systems allow the
simplification of quantum logic gates lanyon2009simplifying , and permit the
enhanced fault tolerance campbell2014enhanced  as well as the efficient
distillation of resource states campbell2012magic  in quantum computation. In
this regard, the capability to sufficiently store the qudit resources with
high efficiency is of crucial importance for constituting high-dimensional
networks so as to distribute high-capacity information in long-distance
quantum communication and facilitate the complex quantum computation.

Qubit memories have been widely demonstrated in many schemes that usually
encode photons in polarization ding2015raman ; vernaz2018highly ;
wang2019efficient ; gundougan2012quantum ; PhysRevLett.111.240503  degree of
freedom (DOF). However, such DOF can only support the two-dimensional
encodings involved with the quantum memory operation. To build up a qudit
memory that can store high-dimensional information, alternative DOFs, such as
which-path chou2005measurement ; moehring2007entanglement ; choi2008mapping ;
pu2017experimental ; PhysRevLett.119.130505 , and time-bin simon2007quantum ;
saglamyurek2011broadband ; clausen2011quantum , have been proposed in a
variety of physical systems. In addition, the photonic transverse spatial
mode, e.g., orbital angular momentum (OAM) mode krenn2017orbital ;
franke2017optical ; ding2013single ; nicolas2014quantum ; parigi2015storage ;
ding2015quantum ; zhang2016experimental ; wang2021efficient ; Ye2022long , has
attracted rapidly growing interest because of its advance of inherent infinite
dimensionality. The storage of these spatial qutrit states with an efficiency
of 20% using the electromagnetically induced transparency (EIT) scheme
ding2014toward  and efficiency of approximately 30% through the off-resonant
Raman protocol zhang2016experimental ; reim2010towards ; reim2011single  have
been reported. However, to date, the maximum available dimensionality of
quantum memory in experiment is limited to d=3 and their efficiencies are far
below the 50% threshold, largely limiting their practical applications in
quantum information processing. The implementation of quantum memories both
having high efficiency and supporting high dimensions is highly desirable but
remains an open challenge.

FIG. 1: Schematic experimental setup. The qudit signal, encoded in POV mode
via SLM 1 and lens L1, is mapped into the atomic ensemble for subsequent
storage. Here, the signal and control fields are both circularly polarized
($\sigma^{+}$), and the control field is beam expanded to have a waist of 4 mm
to completely cover the signal field at the centre of medium.

There are two main challenges to realizing efficient high-dimensional quantum
memories. The first is to establish a uniform light-matter interface to
achieve identical efficiencies for different spatial modes. The imbalanced
storage efficiencies in storing different spatial modes will significantly
degrade the storage fidelity of the qudit state with the increase of
dimensionality. Taking the experiment using Laguerre-Gaussian (LG) mode as a
case in point, the rapidly scaling of the mode waist in $\sqrt{m}$ ($m$ is the
number of modes) ding2015quantum  will lead to significant differences in
light-matter interactions for different modes, thus largely limiting its
applicability in higher-dimensional quantum storage. The second challenge is
to constitute a highly efficient storage medium capable of storing multiple
modes as many as possible grodecka2012high . To achieve this, one needs to
take into account several physical parameters simultaneously in the storage
process, including the transverse spatial extent of the storage medium, the
waist size of the input modes, and the optical depth (OD) of the medium
wang2019efficient ; cao2020efficient . Therefore, the uniform and efficient
storage of a large number of modes is technically challenging.

Here, we demonstrate a high-dimensional quantum memory working up to a
25-dimensional Hilbert space with a storage efficiency of close to 60%, using
the EIT protocol chaneliere2005storage ; eisaman2005electromagnetically ;
zhang2011preparation ; dai2012holographic ; chen2013coherent ; hsiao2018highly
in a laser-cooled atomic ensemble. Through constituting a highly efficient
spatial-mode-independent light-matter interface where photons are encoded in a
unique perfect optical vortex (POV) mode supple  with invariant transverse
size, we are able to store a 25-dimensional qudit by mapping it onto the 25
balanced spatial modes at the centre of the storage medium, and coherently
retrieve these components with identical efficiencies via a control laser. The
demonstrated high-dimensional quantum memory with high efficiency herein is
promising for high-capacity quantum communication and high-dimensional quantum
information processing.

Model and experimental setup. Our memory scheme based on spatial-mode-
independent light-matter interaction is involved with a three-level
$\Lambda$-type atomic system, where the signal field (with a Rabi frequency
$\Omega_{{\rm p}}$) drives the level $\left|{\rm 1}\right\rangle$ to
$\left|{\rm 3}\right\rangle$ and the control field (with a Rabi frequency
$\Omega_{{\rm c}}$) drives the level $\left|{\rm 2}\right\rangle$ to
$\left|{\rm 3}\right\rangle$ (Fig. 1, dashed circle). The dynamical evolution
of the probe field under the slowly-varying envelope approximation can be
described by the Maxwell equation as follows:

$\left[\frac{1}{c}\frac{\partial}{\partial t}+\frac{\partial}{\partial
z}\right]\Omega_{{\rm p}}=i\frac{D_{e{\rm ff}}\Gamma}{2L}\sigma_{31}$ (1)

where $\Gamma$ denotes the decay rate of $\left|{\rm 3}\right\rangle$, $L$ is
the length of medium, and $\sigma_{31}$ represents the atomic coherence
between levels $\left|{\rm 1}\right\rangle$ and $\left|{\rm 3}\right\rangle$.
$D_{e{\rm ff}}\propto N_{{\rm tr}}g_{31}L$ represents the effective OD of an
atomic ensemble, where we define an effective atomic density $N_{{\rm tr}}$
while considering a structured light field interacts with the storage medium
in the transverse orientation. $g_{31}$ represents the photon-atom coupling
coefficient between $\left|{\rm 1}\right\rangle$ and $\left|{\rm
3}\right\rangle$. It can be observed from Eq. (1) that $D_{e{\rm ff}}$
significantly affects the performance of storage, and we derive the numerical
relation between the storage efficiency and OD by solving the Maxwell-Bloch
equations supple .

For a spatial multi-mode quantum memory, it is necessary to take into account
the effective light-matter interaction volume for different spatial modes.
Here, we focus on the coupling of the structure field with the storage medium
in the cross section, because the transverse extent of the storage medium is a
crucial parameter in determining the capacity of multi-mode memory
grodecka2012high . We assume the atomic ensemble with a Gaussian distribution
of the density in the radial direction $N_{{\rm
tr}}(r)=N_{0}\exp[-r^{2}/(2\sigma_{r}^{2})]$. $N_{0}$ refers to the mean
atomic density, and $\sigma_{r}$ represents the half width of the atomic
ensemble supple . In this work, we propose a scheme to establish a uniform
light-matter interface for the memory of a variety of modes via interacting
the photons encoded in POV mode with the storage medium. Theoretically, such
spatial modes feature identical transverse sizes for different $m$, and thus
they are subject to the same $N_{{\rm tr}}(r)$ of atoms when they undergo the
storage process. The interaction strength between the desired POV modes and
medium is uniform, which manifests as the same $D_{e{\rm ff}}$ and ultimately
contributes to the same storage efficiency for different $m$. Based on this
mechanism, we constitute a spatial-mode-independent quantum memory for the
further implementation of storage of high-dimensional quantum states.

FIG. 2: Performance of spatial multi-mode quantum storage. (a) Measured
absorption spectra for various spatial modes versus the signal detuning from
the atomic resonance $\left|1\right\rangle\to\left|3\right\rangle$, where the
relative computer-controlled holograms loaded on the surface of SLM1 are
illustrated in the top right. (b) Transverse intensity distributions of
various modes recorded at the imaging plane of the second 4-f imaging system
before (left) and after (right) storage. (c) Temporal waveforms of input
(blue) and retrieved (red) pulses with temporal lengths of about 500 ns for
different modes. (d) Storage efficiencies versus the quanta of POV mode. The
shaded area represents the maximum fitted value that has been expected, with a
span of 1 sigma. (e) 25$\times$25 input-retrieved cross-talk matrix formed by
the basis set from $\ell=-12$ to 12.

The experimental set-up for a high-dimensional quantum memory is schematically
depicted in Fig. 1. The qudits encoded in each spatial mode are formed on the
basis of the POV eigenstates $\left|\ell\right\rangle$ ($\ell$ is chosen from
-12 to 12), which is accomplished by means of a Fourier transformation of the
Bessel-Gaussian (B-G) state. In this regard, we initially prepare the B-G
states by projecting the attenuated coherent states at the single-photon level
onto a phase-only spatial light modulator (SLM1) to shape the wave-fronts of
photons (Fig. 1, top left). The phase patterns loaded on the SLM are
programmed by a combination of Bessel and Gaussian functions. Lens L1 acting
as a Fourier transformer is then used to transform the B-G states to the POV
states, which are subsequently mapped into the centre of the atomic medium for
storage with the assistance of a carefully aligned 4-f imaging system.

We next store and retrieve the POV states via the EIT storage protocol in a
rubidium medium. To ensure a high storage efficiency of quantum memory, it is
essential to prepare an optically thick atomic ensemble with a large OD, which
is implemented by using a two-dimensional dark-line magneto-optical trap (MOT)
technique in our work. After a programmable storage time, the signal photons
are retrieved from the memory and sent into a qudit state analyser, including
the other 4-f imaging system consisting of lenses L4 and L5, a Fourier lens
L6, as well as a spatial-mode projector based on SLM2, a single-mode fiber
(SMF) and a single-photon counting module (SPCM), to fully characterize the
output states; see the right panel of Fig. 1.

FIG. 3: Characteristics of high-dimensional storage. (a) Distributions of
storage mode bandwidth for different radial wave vectors $k_{r}=1,5,10$, where
the $k_{r}$ of 5 is used in the context. (b) Qudit states with $d=2$, 5, 10,
15, 20 and 25 (see particular expressions in Ref supple ) versus storage
efficiency. (c) Numerical simulation of two-dimensional fidelity as a function
of storage-efficiency-uniformity $\kappa_{1}$. (d) Theoretical analysis of
fidelity versus $\kappa_{1}$ and $\kappa_{2}$ in the case of qudit with $d=3$.

Performance of multi-mode quantum memory. The key to achieving multi-mode
storage in our scheme is to exploit the mode-independent light-matter
interaction. To confirm the accomplishment of this particular photon-atom
interface, we first measure the absorption spectra for a variety of spatial
modes, i.e. $\ell\in\left\\{-12,-6,0,6,12\right\\}$ by scanning the detuning
of signal from $-2{\rm\pi}\times 30$ to $+2{\rm\pi}\times 30$ MHz, as depicted
in Fig. 2(a). The nearly identical OD ($\sim$200) for various
$\left|\ell\right\rangle$ indicates that the interactions between POV photons
and atoms have hardly any correlation with their mode number, thus allowing
our memory to be capable of carrying multiple spatial modes simultaneously. As
shown in Fig. 2(b), the spatial profiles of POV eigenstates with a mean photon
number of n=0.5 in the transverse orientation are detected by an ICCD camera
(iStar 334T series, Andor) working at the single-photon level. The calculated
high values S of similarity ding2013single  between input and retrieved states
are 99.65%, 99.63%, 99.65%, 99.61% and 99.54% for $\ell=-12,-6,0,6,12$
respectively, implying a faithful quantum storage for POV states.

Figure 2(c) shows the temporal waveforms of the input (blue) and retrieved
pulses (red) after a one-pulse-delay storage time for various spatial modes.
As can be seen, the retrievals have almost the same waveforms for different
inputs, providing a clear evidence that our memory exhibits identical
characteristics for different POV modes. To fully analyze the capacity of this
spatial multi-mode quantum memory, we investigate the memory efficiencies of
POV eigenstates across the entire range (from -12 to 12) with a step of
$\Delta\ell=1$; see Fig. 2(d). Their approximately the same values at around
57% clearly illustrate that our memory enables 25 spatial-mode storage with
efficiency beyond 50%. Note that the overall storage-efficiency distributions
for different radial wave vector $k_{r}$ supple  in a wider mode range are
shown in Fig. 3(a). Figure 2(e) gives the experimental cross-talk between the
25 orthogonal bases after retrieval. The average contrast supple , given by
$C=1/25\sum\nolimits_{m}C_{m}$, is estimated to be 92.4$\pm$1.6%, thereby
revealing a low overlap noise between orthogonal spatial modes.

FIG. 4: Demonstration of the storage of quantum states programmed by arbitrary
quanta. (a) Reconstructed real and imaginary parts of density matrices of the
retrieved arbitrary quantum states with $d$=2 in different subspaces. (b)
Single-photon interference fringes for different states. (c) Reconstructed
density matrices of the retrieved qudits with $d$=3 in arbitrarily selected
subspaces. The upper panel illustrates the spatial profiles of the
corresponding quantum states before and after storage. The mean number of
photons per pulse here is $n=0.5$.

In multi-mode memory, the uniform storage efficiency for each POV eigenstate
plays a crucial role in high-dimensional storage. We consider a high-
dimensional quantum superposition state with the dimensionality of $d$, i.e.
the so-called qudit state $\left|{\rm\psi}\right\rangle_{{\rm
Input}}=1/\sqrt{d}\left(\left|\ell_{1}\right\rangle+\left|\ell_{2}\right\rangle+\cdots+\left|\ell_{d}\right\rangle\right)$
as input. The retrieved state after storage can be written as

$\left|\psi\right\rangle_{{\rm
Retrieval}}=1/\sqrt{\sum\nolimits_{m=1}^{d}\eta_{m}^{2}}(\eta_{1}\left|\ell_{1}\right\rangle+\eta_{2}\left|\ell_{2}\right\rangle+\cdots$

$+\eta_{d}\left|\ell_{d}\right\rangle)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ $ (2)

where $\eta_{1},\cdots,\eta_{d}$ denote the storage efficiency for the
corresponding eigenmodes. $\left|\psi\right\rangle_{{\rm Retrieval}}$ can be
further simplified to
$\eta/\sqrt{d}\left(\left|\ell_{1}\right\rangle+\left|\ell_{2}\right\rangle+\cdots+\left|\ell_{d}\right\rangle\right)$
if $\eta_{1},\cdots,\eta_{d}$ are all equal to a constant represented by
$\eta$. In this case, the storage efficiency of the qudits has no dependence
on the dimensionality $d$, as displayed by the results in Fig. 3(b). Thus, our
memory allows storing arbitrarily dimensional qudits with the same efficiency
even when $d$ is up to 25.

Storage fidelity is a critical performance parameter that has to be taken into
account in quantum memory. For the storage of qudit in terms of multiple
spatial modes, its fidelity is extremely sensitive to the uniformity of the
storage efficiency for the internal orthogonal states. For simplicity, we
consider the case of quantum states with $d$ = 2, as shown in Fig. 3(c), where
a parameter $\kappa_{1}$ is defined as the ratio of storage efficiencies
between $\left|\ell_{1}\right\rangle$ and $\left|\ell_{2}\right\rangle$, i.e.
$\kappa_{1}=\eta_{2}/\eta_{1}$. It can be found that the imbalanced atomic
storage ($\kappa_{1}\ll 1)$ would largely reduce the fidelity, as estimated by
the formula $F=\left[{\rm Tr}\left(\sqrt{\sqrt{\rho_{{\rm T}}}\rho_{{\rm
retrieval}}\sqrt{\rho_{{\rm T}}}}\right)\right]^{2}$, where $\rho_{{\rm T}}$
and $\rho_{{\rm retrieval}}$ represent the density matrices corresponding to
the target and retrieval states. In Fig. 4(a), we reconstruct the retrieved
density matrices using the quantum state tomography (QST) method for a set of
qubit states constituted by arbitrary eigenstates (e.g.
$\left|0\right\rangle$, $\left|12\right\rangle$, $\left|5\right\rangle$,
$\left|6\right\rangle$ are chosen herein) after storage. The average fidelity
of 95.8% without any corrections is in good agreement with the theoretical
expectation, and the measured single-photon interference fringes [Fig. 4(b)]
with an average visibility of 92.3% demonstrate that the coherence between two
components of the qubits is well preserved during storage.

In analogy to the case of $d$ = 2, Fig. 3(d) illustrates the effect of
efficiency-uniformity between internal modes on the fidelity for $d$ = 3,
where $\kappa_{2}$ is defined as $\eta_{3}/\eta_{1}$. To obtain a high
fidelity, $\kappa_{1}$ and $\kappa_{2}$ should both approach unity. In Fig.
4(c), we randomly choose three eigenvectors in the range from
$\left|-12\right\rangle$ to $\left|12\right\rangle$ to prepare the high-
dimensional states for storage. The high mean fidelity is measured to be 96.4%
owing to $\kappa_{1}\approx\kappa_{2}\approx 1$. Note that these results can
hardly be obtained in those experiments ding2015quantum  using conventional
vortex modes (e.g., LG mode) because of the inevitable non-uniform efficiency
for different spatial modes. Moreover, we characterize the retrieved state of
$\left|{\rm\psi}_{2}\right\rangle$ for $d$ =5, and the raw fidelity reaches
$90.7\pm 0.7\%$ (the error bar is estimated from Poissonian statistics and
using Monte Carlo simulations), as shown in the right panel of Fig. 5. All
these experimental results indicate our memory capability of storing
arbitrary-mode-encoded qudit states programmed from 25 eigenvectors.

To further prove the quantum nature of the memory, we compare the fidelities
obtained in our experiment with the maximum available fidelities in a
classical memory device based on a completely classical strategy
gundougan2012quantum ; nicolas2014quantum ; parigi2015storage ;
vernaz2018highly . After considering the Poissonian statistics of photon
number for a coherent state, the classical fidelity threshold for a state with
a fixed photon number N can be written as

$F_{{\rm
class}}(n)=\sum\limits_{N=1}^{\infty}\left(\frac{N+1}{N+2}\right)\frac{e^{-n}n^{N}}{(1-e^{-n})N!}$
(3)

where $n$ is the mean photon number per pulse. As presented in Fig. 5, the
solid line is the theoretically classical limit after taking $\eta=0.57$ in
our work. We observe that all the experimental points exceed the classical
benchmark for different mean photon numbers, which confirms the quantum
character of our device.

FIG. 5: Storage in high-dimensional space exceeding the classical benchmark.
The measured fidelities as a function of the mean photon number per pulse n.
The purple/yellow points are experimental data without/with background
subtraction. The blue solid line is the classical limit after considering the
finite storage efficiency and Poissonian statistics of the input.

We now turn to study the capability of our memory to store a 25-dimensional
quantum state. The main challenge to achieving the storage of a 25-dimensional
qudit state is to preserve the identical memory efficiency for each mode, thus
preventing the decay of coherence between 25 spatial modes during the storage
process. Here, a 25-dimensional qudit state $\left|\Psi\right\rangle$ given by
a coherent superposition of 25 individual spatial modes from
$\left|{\rm-12}\right\rangle$ to $\left|{\rm 12}\right\rangle$ is prepared for
the demonstration of 25-dimensional qudit storage, which is represented as

$\left|\Psi\right\rangle=\frac{1}{\sqrt{25}}\sum\limits_{\ell=-12}^{+12}\left|\ell\right\rangle$
(4)

To fully characterize the retrieved state, we perform the high-dimensional QST
supple ; thew2002qudit , where the real and imaginary parts of the
reconstructed density matrix without (with) background correction are plotted
in the logical basis of
{$\left|{\rm-12}\right\rangle,$$\left|{\rm-11}\right\rangle$,
$\left|{\rm-10}\right\rangle$, $\cdots$, $\left|{\rm 12}\right\rangle$}, as
shown in Fig. 6(a,b) (Fig. 6(c,d)), respectively. The raw fidelity between the
retrieved states and ideal state is estimated to be $72.8\pm 0.6\%$, where the
imperfection fidelity is mainly caused by the dark counts of the detector and
residual control laser leakage. After the subtraction of the background, the
fidelity reaches $90.3\pm 0.6\%$, far exceeding the classical limit of 70.2%
for mean photon number $n=0.5$, where the memory efficiency of state
$\left|{\rm\psi}_{6}\right\rangle$ equals to 60% is taken into account. Note
that the residual fidelity is primarily due to the imperfections in the qudit
preparation and measurement. All the above results clearly beat the classical
benchmark, thus demonstrating the quantum character of our 25-dimensional
memory implementation.

FIG. 6: Experimental realization of 25-dimensional qudit storage. The
characterization of the retrieved qudit state
$\left|{\rm\psi}_{6}\right\rangle$ after the storage process by performing
QST. (a)/(c) and (b)/(d) are the real and imaginary parts of the reconstructed
density matrices for retrieved state $\left|{\rm\psi}_{6}\right\rangle$
without/with background correction, respectively.

Conclusion. In summary, we have experimentally demonstrated the efficient
quantum storage for high-dimensional quantum states with d up to 25 using the
POV modes of photons. The reported high-dimensional quantum memory achieves a
storage efficiency of >50%, exceeding the threshold value for practical
quantum information applications. Remarkably, the dimensionality of this
memory is scalable to as high as 100 through further optimization of the waist
of POV modes supple , thus presenting a clear route to the scalability of
dimensions. In addition, our multi-mode memory is also promising for the
compatibility with fiber-based quantum information transfer systems, which are
capable of spatially-structured photon transmission liu2020multidimensional ;
cao2020distribution . The high-dimensional quantum memory demonstrated herein
gives a great perspective for the practical high-capacity and long-distance
quantum communication networks.

This work was supported by National Key R&D Program of China (Grants No.
2017YFA0304800), Anhui Initiative in Quantum Information Technologies (Grant
No. AHY020200), the National Natural Science Foundation of China (Grants No.
U20A20218, No. 61722510, No. 11934013, No. 11604322, No. 12204461), and the
Innovation Fund from CAS, and the Youth Innovation Promotion Association of
CAS under Grant No. 2018490.

## References

  * (1) A. I. Lvovsky, B. C. Sanders, W. Tittel, Optical quantum memory. Nat. photon. 3, 706–714 (2009).
  * (2) N. Sangouard, C. Simon, H. De Riedmatten, N. Gisin, Quantum repeaters based on atomic ensembles and linear optics. Rev. Mod. Phys. 83, 33 (2011).
  * (3) H.-J. Briegel, W. Dür, J. I. Cirac, P. Zoller, Quantum repeaters: the role of imperfect local operations in quantum communication. Phys. Rev. Lett. 81, 5932 (1998).
  * (4) H. J. Kimble, The quantum internet. Nature 453, 1023–1030 (2008).
  * (5) P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, G. J. Milburn, Linear optical quantum computing with photonic qubits. Rev. Mod. Phys. 79, 135 (2007).
  * (6) P. Vernaz-Gris, K. Huang, M. Cao, A. S. Sheremet, J. Laurat, Highly-efficient quantum memory for polarization qubits in a spatially-multiplexed cold atomic ensemble. Nat. Commun. 9, 363 (2018).
  * (7) Y. Wang, et al., Efficient quantum memory for single-photon polarization qubits. Nat. Photon. 13, 346–351 (2019).
  * (8) J. Guo, et al., High-performance raman quantum memory with optimal control in room temperature atoms. Nat. Commun. 10, 1–6 (2019).
  * (9) F. Grosshans, P. Grangier, Quantum cloning and teleportation criteria for continuous quantum variables. Phys. Rev. A 64, 010301(R) (2001).
  * (10) M. Varnava, D. E. Browne, T. Rudolph, Loss tolerance in one-way quantum computation via counterfactual error correction. Phys. Rev. Lett. 97, 120501 (2006).
  * (11) M. Erhard, R. Fickler, M. Krenn, A. Zeilinger, Twisted photons: new quantum perspectives in high dimensions. Light Sci. Appl. 7, 17146–17146 (2018).
  * (12) M. Erhard, M. Krenn, A. Zeilinger, Advances in high-dimensional quantum entanglement. Nat. Rev. Phys. 2, 365–381 (2020).
  * (13) M. Krenn, et al., Generation and confirmation of a (100$\times$100)-dimensional entangled quantum system. PNAS 111, 6243–6247 (2014).
  * (14) H. Bechmann-Pasquinucci, W. Tittel, Quantum cryptography using larger alphabets. Phys. Rev. A 61, 062308 (2000).
  * (15) H. Bechmann-Pasquinucci, A. Peres, Quantum cryptography with 3-state systems. Phys. Rev. Lett. 85, 3313 (2000).
  * (16) N. J. Cerf, M. Bourennane, A. Karlsson, N. Gisin, Security of quantum key distribution using d-level systems. Phys. Rev. Lett. 88, 127902 (2002).
  * (17) S. Walborn, D. Lemelle, M. Almeida, P. S. Ribeiro, Quantum key distribution with higher-order alphabets using spatially encoded qudits. Phys. Rev. Lett. 96, 090501 (2006).
  * (18) L. Sheridan, V. Scarani, Security proof for quantum key distribution using qudit systems. Phys. Rev. A 82, 030301(R) (2010).
  * (19) S. Ecker, et al., Overcoming noise in entanglement distribution. Phys. Rev. X 9, 041042 (2019).
  * (20) B. P. Lanyon, et al., Simplifying quantum logic using higher-dimensional hilbert spaces. Nat. Phys. 5, 134–140 (2009).
  * (21) E. T. Campbell, Enhanced fault-tolerant quantum computing in d-level systems. Phys. Rev. Lett. 113, 230501 (2014).
  * (22) E. T. Campbell, H. Anwar, D. E. Browne, Magic-state distillation in all prime dimensions using quantum reed-muller codes. Phys. Rev. X 2, 041021 (2012).
  * (23) D.-S. Ding, et al., Raman quantum memory of photonic polarized entanglement. Nat. Photon. 9, 332–338 (2015).
  * (24) M. Gündoğan, P. M. Ledingham, A. Almasi, M. Cristiani, H. De Riedmatten, Quantum storage of a photonic polarization qubit in a solid. Phys. Rev. Lett. 108, 190504 (2012).
  * (25) Z. Xu, et al., Long lifetime and high-fidelity quantum memory of photonic polarization qubit by lifting zeeman degeneracy. Phys. Rev. Lett. 111, 240503 (2013).
  * (26) C.-W. Chou, et al., Measurement-induced entanglement for excitation stored in remote atomic ensembles. Nature 438, 828–832 (2005).
  * (27) D. L. Moehring, et al., Entanglement of single-atom quantum bits at a distance. Nature 449, 68–71 (2007).
  * (28) K. S. Choi, H. Deng, J. Laurat, H. Kimble, Mapping photonic entanglement into and out of a quantum memory. Nature 452, 67–71 (2008).
  * (29) Y. Pu, et al., Experimental realization of a multiplexed quantum memory with 225 individually accessible memory cells. Nat. Commun. 8, 1–6 (2017).
  * (30) L. Tian, Z. Xu, L. Chen, W. Ge, H. Yuan, Y. Wen, S. Wang, S. Li, H. Wang, Spatial multiplexing of atom-photon entanglement sources using feedforward control and switching networks. Phys. Rev. Lett. 119, 130505 (2017).
  * (31) C. Simon, H. de. Riedmatten, M. Afzelius, N. Sangouard, H. Zbinden, N. Gisin, Quantum repeaters with photon pair sources and multimode memories. Phys. Rev. Lett. 98, 190503 (2007).
  * (32) E. Saglamyurek, et al., Broadband waveguide quantum memory for entangled photons. Nature 469, 512–515 (2011).
  * (33) C. Clausen, et al., Quantum storage of photonic entanglement in a crystal. Nature 469, 508–511 (2011).
  * (34) M. Krenn, M. Malik, M. Erhard, A. Zeilinger, Orbital angular momentum of photons and the entanglement of laguerre–gaussian modes. Phil. Trans. R. Soc. A. 375, 20150442 (2017).
  * (35) S. Franke-Arnold, Optical angular momentum and atoms. Phil. Trans. R. Soc. A. 375, 20150435 (2017).
  * (36) D.-S. Ding, Z.-Y. Zhou, B.-S. Shi, G.-C. Guo, Single-photon-level quantum image memory based on cold atomic ensembles. Nat. Commun. 4, 2527 (2013).
  * (37) A. Nicolas, et al., A quantum memory for orbital angular momentum photonic qubits. Nat. Photon. 8, 234–238 (2014).
  * (38) V. Parigi, et al., Storage and retrieval of vector beams of light in a multiple-degree-of-freedom quantum memory. Nat. Commun. 6, 7706 (2015).
  * (39) D.-S. Ding, W. Zhang, Z.-Y. Zhou, S. Shi, G.-Y. Xiang, X.-S. Wang, Y.-K. Jiang, B.-S. Shi, G.-C. Guo, Quantum storage of orbital angular momentum entanglement in an atomic ensemble. Phys. Rev. Lett. 114, 050502 (2015).
  * (40) W. Zhang, et al., Experimental realization of entanglement in multiple degrees of freedom between two quantum memories. Nat. Commun. 7, 1–7 (2016).
  * (41) C. Wang, et al., Efficient quantum memory of orbital angular momentum qubits in cold atoms.Quantum Sci. Technol. 6, 045008 (2021).
  * (42) Y.-H. Ye, et al., Long-lived storage of orbital angular momentum quantum states. arXiv preprint arXiv:2203.10884 (2022).
  * (43) D.-S. Ding, et al., Toward high-dimensional-state quantum memory in a cold atomic ensemble. Phys. Rev. A 90, 042301 (2014).
  * (44) K. Reim, et al., Towards high-speed optical quantum memories. Nat. Photon. 4, 218–221 (2010).
  * (45) K. F. Reim, P. Michelberger, K. C. Lee, J. Nunn, N. K. Langford, I. A. Walmsley, Single-photon-level quantum memory at room temperature. Phys. Rev. Lett. 107, 053603 (2011).
  * (46) A. Grodecka-Grad, E. Zeuthen, A. S. Sørensen, High-capacity spatial multimode quantum memories based on atomic ensembles. Phys. Rev. Lett. 109, 133601 (2012).
  * (47) M. T. Cao, et al., Efficient reversible entanglement transfer between light and quantum memories. Optica 7, 1440–1444 (2020).
  * (48) T. Chanelière, et al., Storage and retrieval of single photons transmitted between remote quantum memories. Nature 438, 833–836 (2005).
  * (49) M. D. Eisaman, et al., Electromagnetically induced transparency with tunable single-photon pulses. Nature 438, 837–841 (2005).
  * (50) H. Zhang, et al., Preparation and storage of frequency-uncorrelated entangled photons from cavity-enhanced spontaneous parametric downconversion. Nat. Photon. 5, 628–632 (2011).
  * (51) H.-N. Dai, et al., Holographic storage of biphoton entanglement. Phys. Rev. Lett. 108, 210501 (2012).
  * (52) Y.-H. Chen, M.-J. Lee, I-C. Wang, S. Du, Y.-F. Chen, Y.-C. Chen, I. A. Yu, Coherent optical memory with high storage efficiency and large fractional delay. Phys. Rev. Lett. 110, 083601 (2013).
  * (53) Y.-F. Hsiao, et al., Highly efficient coherent optical memory based on electromagnetically induced transparency. Phys. Rev. Lett. 120, 183602 (2018).
  * (54) See Supplemental Material for the details.
  * (55) R. T. Thew, K. Nemoto, A. G. White, W. J. Munro, Qudit quantum-state tomography. Phys. Rev. A 66, 012303 (2002).
  * (56) J. Liu, et al., Multidimensional entanglement transport through single-mode fiber. Sci. Adv. 6, eaay0837 (2020).
  * (57) H. Cao, et al., Distribution of high-dimensional orbital angular momentum entanglement over a 1 km few-mode fiber. Optica 7, 232–237 (2020).


Neural networks are rapidly gaining interest in nonlinear system identification due to the model's ability to capture complex input-output relations directly from data.
However, despite the flexibility of the approach, there are still concerns about the safety of these models in this context, as well as the need for large amounts of potentially expensive data. Aluminum electrolysis is a highly nonlinear production process, and most of the data must be sampled manually, making the sampling process expensive and infrequent. In the case of infrequent measurements of state variables, the accuracy and open-loop stability of the long-term predictions become highly important. Standard neural networks struggle to provide stable long-term predictions with limited training data. In this work, we investigate the effect of combining concatenated skip-connections and the sparsity-promoting $\ell_1$ regularization on the open-loop stability and accuracy of forecasts with short, medium, and long prediction horizons. The case study is conducted on a high-dimensional and nonlinear simulator representing an aluminum electrolysis cell's mass and energy balance. The proposed model structure contains concatenated skip connections from the input layer and all intermittent layers to the output layer, referred to as InputSkip. $\ell_1$ regularized InputSkip is called sparse InputSkip. The results show that sparse InputSkip outperforms dense and sparse standard feedforward neural networks and dense InputSkip regarding open-loop stability and long-term predictive accuracy. The results are significant when models are trained on datasets of all sizes (small, medium, and large training sets) and for all prediction horizons (short, medium, and long prediction horizons.)

§ INTRODUCTION

There is increasing interest in using machine learning-based methods to develop predictive models directly from data.
Compared to standard system identification methods, the advantage of data-driven modeling is that it does not require any assumptions about the system. However, all phenomena well represented by the data can often be captured accurately.
One example of such a method that has seen widespread popularity in recent years is the NN, which is known to be a universal function approximator.
These are often used in RL to represent a value function or a model for some dynamical system.
However, this approach requires many data points to train effective models, which can be expensive in many domains.
One hypothesis is that NNs are typically overparameterized and require many steps to adjust all parameters.
However, overtraining on the same limited dataset will cause the model to overfit the training data and perform poorly on unseen data.
While overparameterization has been found to aid convergence during training <cit.>, it also introduces redundant information into the weights.

Recent research has found that sparser networks may be the key to training models that generalize across many situations.
In particular, it has been shown empirically that for any dense architecture, there is a high probability that there is a sparse subnetwork that will train faster and generalize better than the full model <cit.>.
This phenomenon is known as the Lottery Ticket Hypothesis. Many sparsification methods can be seen as attempts to extract such a “winning lottery ticket" from an initially dense network.
In system identification, previous work shows that sparsity-promoting $\ell_1$ regularization can benefit model generalization, interpretability, and stability <cit.>.
Group sparsity methods have also been applied to Bayesian recurrent neural networks, with favorable results <cit.>.
There have been numerous advances in this field, and we refer to <cit.> for a recent and comprehensive review.
This work uses the well-known $\ell_1$ regularization to induce sparsity in neural networks.

Another challenge related to using NNs is the choice of architecture and hyperparameters.
Typical networks have multiple layers which are densely connected, although this can vary between domains.
Choosing an appropriate architecture is an art involving trial and error to improve performance and avoid overfitting.
It is commonly understood that the early layers of a neural network significantly impact the overall performance of the network. However, deep networks often suffer from the vanishing or exploding gradient problem, which prevents effective training of these early parameters <cit.>.
Skip-connections were originally proposed to circumvent this by introducing a shorter path between the early layers and the output <cit.>.
They were found to enable the training of significantly deeper networks but may also improve training convergence <cit.>.

In the dynamical systems and control field, models are often designed with a purpose in mind, such as designing a control system or state observer.
Crucially, we are interested in the behavior and performance of the controlled system regarding objectives such as energy efficiency or yield, implying that the model does not need to be perfectly accurate for the entire state space so long as the resulting closed-loop performance is sufficient (known as identification for control (I4C)).
If high-frequency measurements from the system are available, only the short-term behavior of the model is important since any drift out of the operational space is quickly corrected.
However, if measurements are rarely available, such as in the aluminum electrolysis process that we consider, the long-term model behavior and open-loop stability become much more critical.
Stable long-term predictions can be important for decision-making, meaning a model with good long-term stability and accuracy is inherently meaningful.

In this work, we investigate the effects of adding skip connections and $\ell_1$ regularization on the accuracy and stability of these models for short, medium, and long horizons. The following questions are addressed:

* How do skip connections affect the stability and generalization error of neural networks trained on high-dimensional nonlinear dynamical systems?

* How does sparsity affect stability and generalization error for neural networks with skip connections when modeling nonlinear dynamics?

* How does the amount of training data affect neural networks with skip connections compared to neural networks without skip connections?

We make the following contributions:

* We perform a black box system identification of an aluminum electrolysis cell using different NN architectures.

* We demonstrate that the accuracy and open-loop stability of the resulting models is greatly improved by using $\ell_1$ weight regularization and incorporating skip connections into the architecture.

* This advantage is consistent across datasets of varying sizes.

§ THEORY

§.§ Physics-based model for aluminum extraction

Schematic of the setup

NNs are first trained on synthetic data generated from a known PBM.
The model used in this work describes the internal dynamics of an aluminum electrolysis cell based on the Hall-Héroult process.
Fig. <ref> shows a diagram of the electrolysis cell.
Traditional PBMs of such systems are generally constructed by studying the mass/energy balance of the chemical reactions.
The system is described by a set of ODE:

\begin{equation}
\Dot{\state} = \source(\state, \xinput),
\label{eq:nonlin_state_space}
\end{equation}

where $\state \in \mathbb{R}^8$ and $\xinput \in \mathbb{R}^5$ represent the time-varying states and inputs of the system respectively.
The full set of equations are:

\begin{align}
% \begin{aligned}
\Dot{x}_1 &=\,\frac{k_1(g_1 - x_7)}{x_1 k_0} - k_2 (x_6 - g_1) \\
\Dot{x}_2 &=\, u_1 - k_3 u_2\\
\Dot{x}_3 &=\, u_3 - k_4 u_1\\
\Dot{x}_4 &=\, -\frac{k_1 (g_1 - x_7)}{x_1 k_0} + k_2 (x_6 - g_1) + k_5 u_1\\
\Dot{x}_5 &=\, k_6 u_2 - u_4 \\
\Dot{x}_6 &=\, \frac{\alpha}{x_2+x_3+x_4} \Bigg[ u_2 g_5 + \frac{u_2^2 u_5}{2620 g_2} - k_7 (x_6 - g_1)^2 \\
& \, \;  + k_8 \frac{(x_6 - g_1)(g_1 - x_7)}{k_0x_1} - k_9 \frac{x_6 - x_7}{k_{10} + k_{11} k_0  x_1} \Bigg] \nonumber\\
\Dot{x}_7 &=\, \frac{\beta}{x_1} \biggl[ \frac{k_{9} (g_1 - x_7)}{k_{15}k_0 x_1} - k_{12}(x_6 - g_1)(g_1 - x_7)  \\
& \, \; +  \frac{k_{13}(g_1 - x_7)^2}{k_0x_1}
\,  - \frac{x_7 - x_8}{k_{14} + k_{15}k_0  x_1} \biggr]\nonumber\\
\Dot{x}_8 &=\, k_{17} k_9 \left(\frac{x_7 - x_8}{k_{14} + k_{15} k_0 \cdot x_1} - \frac{x_8 - k_{16}}{k_{14} + k_{18}}\right),
% \end{aligned}
\end{align}

where the intrinsic properties $g_i$ of the bath mixture are given as:

\begin{align}
g_1 &=  991.2 + 112 c_{x_3} + 61 c_{x_3}^{1.5} - 3265.5 c_{x_3}^{2.2} \label{eq:liquidus_temp}  \\
&- \frac{793 c_{x_2}}{- 23 c_{x_2} c_{x_3} - 17 c_{x_3}^{2} + 9.36 c_{x_3} + 1} \nonumber \\
g_2 &=  \text{exp}\,\left(2.496 - \frac{2068.4}{273+x_6} - 2.07c_{x_2}\right)\\
g_3 &= 0.531 + 3.06 \cdot 10^{-18} u_{1}^{3} - 2.51 \cdot 10^{-12} u_{1}^{2} \\ &+ 6.96 \cdot 10^{-7} u_{1}
- \frac{14.37 (c_{x_2} - c_{x2,crit}) - 0.431}{735.3 (c_{x_2} - c_{x2,crit}) + 1}  \nonumber\\
g_4 &= \frac{0.5517 + 3.8168 \cdot 10^{-6}u_2}{1 + 8.271 \cdot 10^{-6}u_2}\\
g_5 &= \frac{3.8168\cdot 10^{-6} g_3 g_4 u_2}{g_2(1 -g_3)}.
\end{align}

See Table <ref> for a description of these quantities.
The dynamics of the system are relatively slow. The control inputs $u_1, \; u_3$ and $u_4$ are therefore well modeled as impulses representing discrete events involving the addition or removal of substances. This results in step changes in the linear states $x_2, x_3, x_5$, which act as accumulator states for the mass of the corresponding substance (see Table <ref>). The control inputs $u_2$ and $u_5$ are piecewise constant and nonzero. The inputs $\xinput$ are determined by a simple proportional controller $\boldsymbol{\pi}(\state)$. The simulation model is derived in <cit.>, and we refer to that article for the values of the simulation parameters and further details.

Table of states, inputs, and other quantities used to model the electrolysis cell

Variable     Physical meaning     Units
$x_1$     Mass side ledge     $\si{kg}$

$x_2$     Mass Al2O3     $\si{kg}$

$x_3$     Mass AlF_3     $\si{kg}$

$x_4$     Mass Na_3 AlF_6     $\si{kg}$

$x_5$     Mass metal     $\si{kg}$

$x_6$     Temperature bath     $\si{\degreeCelsius}$

$x_7$     Temperature side ledge     $\si{\degreeCelsius}$

$x_8$     Temperature side wall     $\si{\degreeCelsius}$

$u_1$      Al2O3 feed    $\si{kg/s}$

$u_2$     Line current      $\si{kA}$

$u_3$     AlF_3 feed     $\si{kg/s}$

$u_4$     Aluminum tapping     $\si{kg/s}$

$u_5$     Anode-cathode distance     $\si{cm}$

$c_{x_2}$     Al2O3 mass ratio $x_2/(x_2 + x_3 + x_4)$      -

$c_{x_3}$     AlF_3   mass ratio $x_3/(x_2 + x_3 + x_4)$     -

$g_1$     Liquidus temperature        $\si{\degreeCelsius}$

$g_2$     Electrical conductivity     $\si{\siemens\meter}$

$g_3$     Bubble coverage             -

$g_4$     Bubble thickness            $\si{cm}$

$g_5$     Bubble voltage              $\si{\volt}$

§.§ Deep neural network with skip connections

A NN with $L$ layers can be compactly written as an alternating composition of affine transformations $\mathbf{Wz + b}$ and nonlinear activation functions $\boldsymbol{\sigma} : \mathbb{R}^n \mapsto \mathbb{R}^n$:

\begin{equation}
\begin{aligned}
\mathbf{z}_{i} = \boldsymbol{\sigma}_i\left(\mathbf{W}_i\mathbf{z}_{i-1} + \mathbf{b}_i\right)
\end{aligned}
\end{equation}

where $\mathbf{z}_0$ is the input to the network, the activation function $\boldsymbol{\sigma}_i$, weight matrix $\mathbf{W}_i$, and bias vector $\mathbf{b}_i$ correspond to the $i$th layer of the network. The universal approximation property of NNs makes them very attractive as a flexible model class when a lot of data is available. The representation capacity is generally understood to increase with both the depth and the width (the number of neurons in each layer), although early attempts to train very deep networks found them challenging to optimize using backpropagation due to the vanishing gradients problem.
One of the major developments that enabled researchers to train deep NNs with many layers is the skip connection. A skip connection is simply an additional inter-layer connection that bypasses some of the layers of the network. This provides alternate pathways through which the loss can be backpropagated to the early layers of the NN, which helps mitigate the issues of vanishing and exploding gradients, which were major hurdles to training deeper models. In this work, we utilize a modified DenseNet architecture as proposed by <cit.>, where the outputs of earlier layers are concatenated to all the consecutive layers. We simplify the structure such that the model only contains skip connections from the input layer to all consecutive layers. We call this architecture InputSkip, which has reduced complexity compared to DenseNet:

\begin{equation} \label{eq:inputskip}
\begin{aligned}
\mathbf{z}_{1} &= \boldsymbol{\sigma}_1(\mathbf{W}_1\mathbf{z}_{0} + \mathbf{b}_1) \\
\mathbf{z}_{i} &= \boldsymbol{\sigma}_i\left(\mathbf{W}_i\begin{bmatrix}\mathbf{z}_{i-1} \\ \mathbf{z}_{0}\end{bmatrix} + \mathbf{b}_i\right) , i>1
\end{aligned}
\end{equation}

The output of each layer (excl. the first) becomes a sum of a linear and a nonlinear transformation of the initial input $\mathbf{x}$. Hence, the skip connections from the input layer to consecutive layers facilitate the reuse of the input features for modeling different linear and nonlinear relationships more independently. The InputSkip architecture with four hidden layers is illustrated in Figure <ref>.

InputSkip architecture with 4 hidden layers

§ METHOD AND SETUP

In this section, we present all the details of data generation, its preprocessing, and the methods required to reproduce the work. The steps can be briefly summarized as follows:

* Use (<ref>) with random initial conditions to generate 140 trajectories with 5000 timesteps each. Set aside 40 for training and 100 for testing. Construct three datasets by selecting 10,20, and 40 trajectories, respectively.

* For each model class and dataset, train ten instances on the training data.

* Repeat all experiments with $\ell_1$ regularization, see loss function in (<ref>).

* Use trained models to generate predicted trajectories along the test set and compare them to the 100 test trajectories.

§.§ Data generation

The state trajectories used in the test and training sets were generated by integrating Equation (<ref>) with the numerical RK4 integration scheme with a fixed timestep $\timestep = \SI{10}{\second}$ on the interval $[0, 5000\timestep]$.
The initial conditions were sampled uniformly from the intervals shown in Table <ref> to generate 140 unique trajectories.
A total of 40 trajectories were set aside for training and 100 of the trajectories as a test set.
The 40 training trajectories were used to create three datasets of varying sizes (small, medium, large), namely 10, 20, and 40 trajectories, containing 50000, 100000, and 200000 individual data points.

Equation (<ref>) also depends on the input signal $\xinput$.
In practice, this is given by a deterministic control policy $\xinput = \boldsymbol{\pi}(\state)$ that stabilizes the system and keeps the state $\state$ within some region of the state space that is suitable for safe operation.
We found that this was insufficient to successfully train our models because the controlled trajectories showed minimal variation after some time, despite having different initial conditions.
This lack of diversity in the dataset resulted in models that could not generalize to unseen states, which frequently arose during evaluation.

To inject more variety into the data and sample states $\state$ outside of the standard operational area, we used a stochastic controller

$$\boldsymbol{\pi}_s(\state) = \boldsymbol{\pi}(\state) + \mathbf{r}(t)$$

that introduced random perturbations $\mathbf{r}(t)$ to the input.
These perturbations were sampled using the APRBS method proposed by <cit.>.

In system identification, it is typical to optimize the model to estimate the function $\mathbf{\dot{x} = f(x,u)}$. However, this is not feasible for (<ref>) because the inputs $\xinput$ are not differentiable. Instead, the trajectories are discretized using the forward Euler difference:

\begin{equation}\label{eq:forward_euler}
\mathbf{y}_k = \frac{\state_{k+1} - \state_{k}}{\timestep}
\end{equation}

The datasets are then constructed as sets of the pairs $([\state_k,\xinput_k], \mathbf{y}_k)$.
In practice, measurements will be noisy and the state trajectories must be estimated using a filtering method, e.g., moving horizon estimation.

Initial conditions intervals for $\state$

Variable      Initial condition interval
$x_1$     $[2060,\;4460]$

$c_{x_2}$     $[0.02,\; 0.05]$

$c_{x_3}$     $[0.09,\; 0.13]$

$x_4$    $[11500,\; 16000]$

$x_5$     $[9550,\; 10600]$

$x_6$     $[940,\; 990]$

$x_7$     $[790,\; 850]$

$x_8$     $[555,\; 610]$

§.§ Model architectures

Two different architectures are evaluated in this case study: a standard feed-forward Multi-Layer Perceptrons (MLP) referred to as PlainNet and the modified MLP with concatenated skip-connections from the input layer called InputSkip, see Fig. <ref> for illustration. Moreover, both structures are trained with and without the sparsity promoting $\ell_1$ regularization, yielding four model structures: PlainDense, PlainSparse, InputSkipDense, and InputSkipSparse. The input layer of each of the models is the concatenation of the measured state $\mathbf{x}_k \in \mathbb{R}^8$, or the estimated state $\hat{\mathbf{x}}_k\in \mathbb{R}^8$ at timestep $k$, and the control input vector $\mathbf{u}_k\in \mathbb{R}^8$ at timestep $k$, yielding a vector $\mathbf{z}_0 = \{\mathbf{x}_k, \; \mathbf{u}_k\} \in  \mathbb{R}^{13}$. Each of the structures has four hidden layers, and each of the layers has 25 neurons. In addition, the input vector $\mathbf{z_0}$ is concatenated to each of the hidden layers in the InputSkip structures, such that each of the hidden layers in the InputSkip structures has $25 + 13 = 38$ states. All models output an estimate of the time derivative of the state variables at timestep $k$ $\Dot{\mathbf{x}}_k \in \mathbb{R}^8$. Each model class's sparse and dense structures start with the same architecture before training, but for sparse structures, many of the neurons and weights are zeroed out by the $\ell_1$ regularization term. Since InputSkip has more states in the hidden layers than Plain structures due to the input vector being concatenated to each layer's output, it is reasonable to ask whether the Plain structures should have more neurons in the hidden layers. This is tested, and it turns out that this does not benefit the structures regarding the evaluation measures.

§.§ Training setup

The models are trained by minimizing the following loss function using stochastic gradient descent:

\begin{equation}
\label{eq:l1}
\mathbf{J}_{\theta} = \frac{1}{|\mathcal{B}|}\sum_{i\in\mathcal{B}} (\mathbf{\mathbf{y}_i}-\mathbf{\net(\state_i, \xinput_i}))^2 + \lambda \sum_{j=1}^{L}|\mathbf{W}_j|
\end{equation}

where the batch $\mathcal{B}$ is a set of indices corresponding to a random subset of examples from the data, $L$ is the number of layers of the NN, and $\lambda$ is the regularization parameter. This loss function is the sum of the MSE of the model $\net$ with respect to the regression variables $\mathbf{y}$, and the $\ell_1$ norm of the connection weight matrices $\mathbf{W}_i$ in all layers. We used a batch size of $|\mathcal{B}| = 128$. We used the popular ADAM solver proposed by <cit.> with default parameters to minimize (<ref>). The dense model structures PlainDense and InputSkipDense were trained with $\lambda=0$, and the sparse model structures PlainSparse and InputSkipSparse were trained with $\lambda=10^{-4}$.

§.§ Evaluation of model accuracy

Starting from a given initial condition $\mathbf{x}(t_0)$, the model $\net(\state, \xinput)$ is used to generate an estimated trajectory using the recurrence:

\begin{equation}
\label{eq:euler_forecast}
\hat{\state}_{k+1} = \hat{\state}_{k} + \timestep\,\net(\hat{\state}_{k}, \xinput_k)
\end{equation}

where $\hat{\state}_{0} = {\state}_{0}$.
Applying multi-step or higher order Runge-Kutta methods to NN models is possible. However, these methods require multiple evaluations of the model per timestep, which increases the computation and memory needed to perform automatic differentiation.
The forward Euler method is preferred, as it only evaluates the model once per timestep.
In the approach outlined here, this does not incur significant discretization errors, as (<ref>) effectively reverses the discretization step in (<ref>).
The input signal $\xinput_k$ is sampled directly from the test trajectory.
Borrowing a term from the field of time-series analysis, this is referred to as a rolling forecast.
To evaluate the accuracy of a model over multiple trajectories, we define the Average Normalized Rolling Forecast Mean Squared Error (AN-RFMSE):

\begin{equation}
\textrm{AN-RFMSE} = \frac{1}{p}\sum_{i=1}^p\frac{1}{n}\sum_{j=1}^n\left(\frac{\hat{x}_i(t_j) - x_i(t_j)}{\std(x_i)}\right)^2,
\label{eq:AN-RFMSE}
\end{equation}

where $\hat{x}_i(t_j)$ is the model estimate of the simulated state variable $x_i$ at time step $t_j$, $\std(x_i)$ is the standard deviation of variable $x_i$ in the training set $\mathcal{S}_{train}$, $p=8$ is the number of state variables and $n$ is the number of time steps being averaged over.

§.§ Evaluation of model stability

A symptom of model instability is that its predictions can blow up, characterized by a rapid (often exponential) increase in prediction error. More precisely, a blow-up is said to occur when all system states' normalized mean absolute error exceeds three (this corresponds to standard deviations):

\begin{equation}
\max_{j<n} \left[\frac{1}{p}\sum_{i=1}^p\left(\frac{\left\vert\hat{x}_i(t_j) - x_i(t_j)\right\vert}{\std(x_i)}\right)\right] > 3
\label{eq:blowup-detection}
\end{equation}

where $p=8$ is again the number of state variables and $n$ is the number of time steps to consider.
Equation (<ref>) is conservative. However, this does not lead to a significant underestimation of the number of blow-ups. This is because once a model starts to drift rapidly, it quickly exceeds the threshold.

§ RESULTS AND DISCUSSIONS

This section reports empirical results for the model accuracy and stability of the different model classes (PlainDense, PlainSparse, InputSkipDense, InputSkipSparse).
A Monte Carlo analysis is performed by training ten instances of each model class and evaluating these on the test set consisting of 100 trajectories, where each trajectory has a length of 5000 timesteps. The test set is generated as described in Section <ref>, using the simulation model in Equation (<ref>). The evaluation procedure follows; the models are given initial values for each state variable trajectory in the test set. Then, the models forecast state values at the consecutive time steps for each test set trajectory as described in (<ref>) without feedback from measurements of the state variables. The resulting forecasts are evaluated according to the prediction accuracy measure in (<ref>) and the forecast stability measure in (<ref>). The prediction accuracy of different model classes is reported in Fig. <ref>, and the empirical model stability results are reported in Fig. <ref>. Accuracy and stability are reported for different forecasting horizons. We repeat the experiments for all model classes trained on three different dataset sizes to study the data efficiency of the models.

Fig. <ref> presents the total number of blow-ups recorded within each model class after $100\timestep$, $2000\timestep$, and $5000\timestep$ (short, medium, and long term respectively).
For simplicity, blow-ups were detected by thresholding the computed variance of a predicted trajectory and manually inspected.
It is clear that for short time horizons, all the models exhibit robust behavior independently of the size of the training datasets. However, for medium and long time horizons, PlainDense, PlainSparse, and InputSkipDense architectures exhibit a significant number of blow-ups and, therefore, instability. Figs. <ref> - <ref> show that PlainDense is generally the most unstable, with up to 41% of all trajectories resulting in a blow-up. For the smallest amount of training data (see Fig. <ref>) PlainSparse and InputSkipDense have similar blow-up frequencies. The PlainSparse architecture shows significantly better stability for larger datasets than both PlainDense and InputSkipDense. InputSkipDense and PlainDense both show better stability with increasing training data regarding fewer blow-ups. However, both these dense models still suffer from high blow-up rates.

In comparison, almost no blow-ups are recorded using the InputSkipSparse architecture, even for the small training dataset. In Figure <ref>, the orange bars corresponding to the blow-up frequency of InputSkipSparse models are not visible for any training sets due to the significantly lower number of blow-ups. For InputSkipSparse models trained on the smallest dataset, only 3 out of 1000 possible blow-ups were reported for the longest horizon. Apart from that, no blow-ups were reported for the InputSkipSparse models.
Only a few blow-ups were recorded after $5000\timestep$ in the medium term.

Trained on smallest dataset with 50000 datapoints

Trained on medium sized dataset with 100000 datapoints

Trained on largest dataset with 200000 datapoints

max width=

[legend columns=4,legend style=draw=none, column sep=1ex,legend entries= PlainDense,PlainSparse,InputSkipDense,InputSkipSparse]
Tpurple!50, fill=Tpurple!50,area legend
Tblue!50, fill=Tblue!50,area legend
Tgreen!50, fill=Tgreen!50,area legend
Torange!50, fill=Torange!50,area legend

Divergence plot: Number of trajectories that blow-up over different time horizons according to the measure defined in Equation (<ref>). The total number of trajectories is 1000, so the values can be read as a permille.

Fig. <ref> presents a violin plot of the accuracy of each model class, expressed in terms of AN-RFMSE over different time horizons.
A larger width of the violin indicates a higher density of that given RFMSE value, while the error bars show the minimum and maximum recorded RFMSE values.
The model estimates that blew up (see Fig. <ref>) are excluded as outliers.
In this way, the generalization performance of the models is estimated only within their regions of stability.
A potential pitfall of excluding these outliers is that model classes that blow up often have their worst scores removed, thus biasing the distribution towards lower scores.
Despite this, low accuracy appears to correlate with a high blow-up rate. The InputSkipSparse architecture is consistently more accurate (up to an order of magnitude) than the others in the long term.

Trained on smallest dataset with 50000 datapoints

Trained on medium sized dataset with 100000 datapoints

Trained on largest dataset with 200000 datapoints

max width=

[legend columns=4,legend style=draw=none, column sep=1ex,legend entries= PlainDense,PlainSparse,InputSkipDense,InputSkipSparse]
Tpurple!50, fill=Tpurple!50,area legend
Tblue!50, fill=Tblue!50,area legend
Tgreen!50, fill=Tgreen!50,area legend
Torange!50, fill=Torange!50,area legend

Model accuracy is expressed in terms of AN-RFMSE over different horizons. AN-RFMSE is defined for a single model forecast of a single trajectory in (<ref>). Ten models of each of the model types (PlainDense, PlainSparse, InputSkipDense, InputSkipSparse) are trained on the 50000 data points in Figure <ref>, 100000 data points in Figure <ref>, and 200000 data points in Figure <ref>. The model estimates that blow up (see Figure <ref>) are excluded. The plot shows that sparse models with skip connections (InputSkipSparse) are consistently more accurate than sparse and dense models without skip connections.

Side ledge mass $x_1$

Alumina mass $x_2$

Aluminum fluoride $x_3$

Molten cryolite $x_4$

Produced aluminum $x_5$

Bath temperature $x_6$

Side ledge temperature $x_7$

Side wall temperature $x_8$

max width=1

[legend columns=6,legend style=draw=none, column sep=2ex,legend entries=Truth,InputSkipSparse,PlainSparse,$99.7\%$ conf. PlainSparse,$99.7\%$ conf. InputSkipSparse]
black,thick, sharp plot
Torange,thick, dashed,sharp plot
Tblue,thick, dashed,sharp plot
Tpurple!30, fill=Tpurple!20, area legend
Torange!40, fill=Torange!40, area legend

Rolling forecast of a representative trajectory from the test set (100 trajectories total)

Fig. <ref> shows the effect of InputSkipSparse compared to PlainSparse for a single representative test-set trajectory and is not meant to the significance of the results. The significance of the results can be found in Fig. <ref> and Fig. <ref>, which show results for the entire test set.

§ CONCLUSION AND FUTURE WORK

This work compared the performance of two different model structures trained with and without sparsity promoting $\ell_1$ regularization. The two model types are standard MLP and a more specialized architecture that includes skip connections from the input layer to all consecutive layers, yielding four different model structures: PlainDense, PlainSparse, InputSkipDense, and InputSkipSparse. The main conclusions of the article are as follows:

* NNs with skip connections are more stable for predictions over long time horizons compared to standard MLPs. Furthermore, the accuracy of NNs with skip connections is consistently higher for all forecasting horizons.

* The application of sparsity-promoting $\ell_1$ regularization significantly improves the stability of the standard MLP and InputSkip architectures. This improvement was more apparent for models with the InputSkip architecture.

* The InputSkipSparse showed satisfactory stability characteristics even when the amount of training data was restricted, suggesting that this architecture is more suitable for system identification tasks than the standard MLP structure.

The case study shows that both sparsity-promoting regularization and skip connections can result in more stable NN models for system identification tasks while requiring fewer data and improving their multi-step generalization for both short, medium, and long prediction horizons.
Despite the encouraging performance of the sparse-skip networks, it is yet to be determined if the benefits also extend to the case of noisy measurements.

This case study also has relevance beyond the current case study.
In more realistic situations, we often have a partial understanding of the system we wish to model (see (<ref>)) and only wish to use data-driven methods to correct a PBM when it disagrees with the observations (e.g., due to a faulty assumption).
As shown in <cit.>, combining PBMs and data-driven methods in this way also has the potential to inject instability into the system.
Finding new ways to improve or guarantee out-of-sample behavior for data-driven methods is therefore paramount to improving such systems' safety.