# Superpolynomial Quantum Enhancement in Polaritonic Neuromorphic Computing

Huawen Xu<EMAIL_ADDRESS>Division of Physics and Applied Physics,
School of Physical and Mathematical Sciences, Nanyang Technological
University, 21 Nanyang Link, Singapore 637371, Singapore  Tanjung Krisnanda
Division of Physics and Applied Physics, School of Physical and Mathematical
Sciences, Nanyang Technological University, 21 Nanyang Link, Singapore 637371,
Singapore  Wouter Verstraelen  Division of Physics and Applied Physics, School
of Physical and Mathematical Sciences, Nanyang Technological University, 21
Nanyang Link, Singapore 637371, Singapore  Timothy C. H. Liew
<EMAIL_ADDRESS>Division of Physics and Applied Physics, School of
Physical and Mathematical Sciences, Nanyang Technological University, 21
Nanyang Link, Singapore 637371, Singapore  MajuLab, International Joint
Research Unit UMI 3654, CNRS, Université Côte d’Azur, Sorbonne Université,
National University of Singapore, Nanyang Technological University, Singapore
Sanjib Ghosh<EMAIL_ADDRESS>Division of Physics and Applied Physics,
School of Physical and Mathematical Sciences, Nanyang Technological
University, 21 Nanyang Link, Singapore 637371, Singapore

###### Abstract

Recent proof-of-principle experiments have demonstrated the implementation of
neuromorphic computing using exciton-polaritons, making use of coherent
classical states [D. Ballarini _et al_., Nano Lett. 20, 3506 (2020)]. At the
same time, it is expected that nonlinear exciton-polaritons can reach a
quantum regime forming non-classical states. Here we consider theoretically
the quantum nature of exciton polaritons and predict a superpolynomial quantum
enhancement in image recognition tasks. This is achieved within experimentally
accessible parameters.

##  I I. INTRODUCTION

The Exciton-polariton has emerged as an attractive quasiparticle, a hybrid of
light (photon) and matter (exciton) in semiconductor microcavities, which
possesses the characteristics of long dephasing time and ultrafast response
inherited from their photonic part as well as strong nonlinearity from their
excitonic part Carusotto and Ciuti (2013); Kasprzak  _et al._ (2006). These
properties have encouraged works developing polaritonic information processing
devices Sanvitto and Kéna-Cohen (2016); Fraser (2017). In particular, the
experimental realization of various processing components has been reported,
including: logic gates Leyder  _et al._ (2007); Cancellieri  _et al._ (2015),
diodes Nguyen  _et al._ (2013), amplifiers Niemietz  _et al._ (2016), switches
Gao  _et al._ (2012); De Giorgi  _et al._ (2012); Dreismann  _et al._ (2016);
Grosso  _et al._ (2014), couplers Rozas  _et al._ (2020); Klaas  _et al._
(2019), junctions Winkler  _et al._ (2017), transistors Ballarini  _et al._
(2013); Lewandowski  _et al._ (2017); Zasedatelev  _et al._ (2019), routers
Flayac and Savenko (2013); Marsault  _et al._ (2015), high speed transmission
lines Lerario  _et al._ (2017) and memories Cerna  _et al._ (2013). Although
individual elements are well studied with exciton-polaritons, very few
experimental works have studied their combination into complete information
processing systems.

Theoretical works have considered that rather than imitating CMOS logic,
exciton-polaritons might be more suited to alternative architectures of
information processing, such as simulators of other systems Kalinin and
Berloff (2019) or Hamiltonian optimizers Kalinin and Berloff (2018), cellular
automata Banerjee and Liew (2020), or neural networks Espinosa-Ortega and Liew
(2015). These proposals have been based on the high degree of tuneability
available in exciton-polariton systems, which is particularly useful in neural
networks where functionality is determined by the strength of weight
connections between network nodes. In general, neural network architectures
serve as a good platform to apply machine learning algorithms, which
circumvent the need to encode the function of the network in its design
directly but rather allow it to be learned. Indeed, simple feedback algorithms
have been used in exciton-polariton experiments to induce specific patterns
Ohadi  _et al._ (2017); Töpfer  _et al._ (2020). However, so far the
implementation of the most well-known feedforward neural networks has not been
attempted with exciton-polaritons, partly due to the large number of weight
connections that would need individual control.

The same challenge of needing to simultaneously control many weight
connections has been avoided in many other hardware realizations of neural
networks by using so-called reservoir computing networks Tanaka  _et al._
(2019). These architectures are based on completely randomly connected
distributions of network nodes Jaeger and Haas (2004); Lukoševičius and Jaeger
(2009); Enel  _et al._ (2016), where all training is implemented in a single
output layer, which also greatly simplifies learning algorithms Lukoševičius
(2012). Reservoir computing was considered theoretically for exciton-polariton
systems Opala  _et al._ (2019); Xu  _et al._ (2020) and implemented
experimentally shortly after Ballarini  _et al._ (2020); Mirek  _et al._
(2021), where benchmarking for the task of recognizing a standard set of hand-
written digits showed the exciton-polariton system to have a higher success
rate than previous hardware implementations of reservoir computing.

All the aforementioned works on exciton-polaritons have operated with
classical coherent states, formed either through direct laser excitation or
the formation of a polariton condensate. On the other hand, polaritons are
quantum particles, which has been evidenced in different experiments. In
particular, the onset of non-classical correlations Savasta  _et al._ (2005)
and squeezing Karr  _et al._ (2004); Boulier  _et al._ (2014) has been
reported, while more recent experiments showed the preservation of polariton
entanglement Cuevas  _et al._ (2018) and imaging of the single polariton
wavefunction Suárez-Forero  _et al._ (2020). Other recent experiments Muñoz-
Matutano  _et al._ (2019); Delteil  _et al._ (2019) have shown the onset of
the polariton blockade Verger  _et al._ (2006), where antibunching is
generated from the interactions between polaritons themselves. Although the
amount of antibunching measured was limited by the weak polariton-polariton
interaction strength, a number of methods exist for increasing it Kyriienko
_et al._ (2020); Rosenberg  _et al._ (2018); Jia  _et al._ (2018); Stefanatos
and Paspalakis (2020).

Here, we consider the performance of polariton reservoir computers accounting
for their quantum nature. In principle, the larger Fock space can be expected
to allow a quantum enhancement, as has been demonstrated for quantum
reservoirs based on spin systems Fujii and Nakajima (2017); Nakajima  _et al._
(2019); Chen and Nurdin (2019); Chen  _et al._ (2020) and a single nonlinear
oscillator Govia  _et al._ (2021) or a qudit Kalfus  _et al._ (2021). As
polariton reservoir computers have previously been benchmarked for image
recognition, we will consider images with their pixels converted into a time-
dependent signal that excites the quantum polariton reservoir with a resonant
amplitude and/or phase modulated laser. The output from the quantum polariton
reservoir is represented by the quasiprobability distribution in phase space,
namely its Wigner function Olivares (2012), which will be used to classify the
input image. By comparing to simulations with classical polariton reservoirs,
we find that the quantum polariton system is able to classify images with a
single polariton node! We will continue to refer to this node as the
“reservoir” as it nevertheless has a significant Fock space and fulfils the
role of the much larger reservoir of classical reservoir computers, which
cannot attain the same error rate with even hundreds of classical nodes. The
single quantum polariton node exceeds also the limit of a classical linear
classifier and is found empirically to have a superpolynomial quantum
enhancement. We also compare possibilities using two-photon excitation, which
allows further enhancement.

##  II II. Quantum polariton reservoir

The physical model considered here consists of only a single polariton mode,
where the Fock states of the quantum mode form the quantum reservoir network,
see Fig. 1a. The bare Hamiltonian $\mathcal{\hat{H}}_{0}$ for the polariton
mode is taken as

$\displaystyle\mathcal{\hat{H}}_{0}=E\hat{a}^{\dagger}\hat{a}+\alpha\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a},$
(1)

where $\hat{a}^{\dagger}$ ($\hat{a}$) is the creation (annihilation) operator,
$E$ is the energy of the polariton mode and $\alpha$ is the Kerr nonlinear
interaction strength. For encoding information, we introduce a coherent single
photon ($\mathcal{SP}$) and two photon ($\mathcal{TP}$) pumping, such that the
total Hamiltonian (in a rotating frame at the frequency of the pump
$\omega_{p}$) is written as

$\displaystyle\mathcal{\hat{H}}$ $\displaystyle=$
$\displaystyle\Delta\hat{a}^{\dagger}\hat{a}+\alpha\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}$
(2)
$\displaystyle+F(\hat{a}^{\dagger}+\hat{a})+P(e^{i\Theta}\hat{a}^{\dagger}\hat{a}^{\dagger}+e^{-i\Theta}\hat{a}\hat{a}),$

where $F$ and $P$ are real amplitudes of the $\mathcal{SP}$ and $\mathcal{TP}$
pumping, respectively, and we allowed for a phase difference of the
$\mathcal{TP}$ pump, denoted by $\Theta$. The detuning term is defined as
$\Delta=E-\hbar\omega_{p}$, and we have assumed that the pumping frequency of
$\mathcal{SP}$ ($\mathcal{TP}$) is $\omega_{p}$ ($2\omega_{p}$). See also Ref.
Bartolo  _et al._ (2016), in which a similar treatment is presented. The
network structure composed of the Fock states can be illustrated as follows.
The single pumping term in Eq. (2) expressed in the Fock basis reads
$\sum_{n=0}F\sqrt{n+1}(|n+1\rangle\langle n|+|n\rangle\langle n+1|)$. This way
the $\mathcal{SP}$ is seen as nearest neighbour connection between the Fock
states. Similarly, the two photon pumping serves as the connection between
next nearest neighbours. In our analysis, we compare four different pumping
scenarios: $\mathcal{SP}$ ($P=0$); $\mathcal{TP}$ ($F,\Theta=0$);
$\mathcal{SPTP}$ ($\Theta=0$); as well as $\mathcal{SPTP}^{\prime}$ where
$\Theta=\Theta(t)\neq 0$.

Figure 1: The schematic for image classification using quantum polariton
reservoir computing. a) Input signals excite a single nonlinear polariton mode
in the form of pumping amplitudes and/or phases. After a certain time $\tau$,
the quasiprobability distribution in the phase space (Wigner function) is
calculated. By applying trained output weights ($\bm{W}^{\text{out}}$) to the
Wigner function, a final output classifies the input signals (images). b) Each
image (from the MNIST standard data set) is first compressed from $28\times
28$ resolution to $4\times 4$. The pixels of the compressed image serve as
time-dependent input signals.

In our simulations, we evolve the density matrix $\rho$ under the different
pumping methods, where the initial state of the polariton mode is given by the
vacuum state. The dynamics of the system can be described by the quantum
master equation:

$\displaystyle
i\hbar{\dot{\rho}}=[\mathcal{\hat{H}},\rho]+\frac{i\gamma}{2}\mathscr{L}(\rho,\hat{a}).$
(3)

where $\mathcal{\hat{H}}$ is the Hamiltonian of Eq. (2),
$\mathscr{L}(\rho,\hat{a})$ is the Lindblad dissipation superoperator defined
as $\mathscr{L}(\rho,\hat{a})\equiv
2\hat{a}\rho\hat{a}^{\dagger}-\hat{a}^{\dagger}\hat{a}\rho-\rho\hat{a}^{\dagger}\hat{a}$,
and $\gamma$ is the dissipation rate. As can be seen from Fig. 1a, the input
information is first fed into the single polariton node. After a finite time
$\tau$ of evolution, the Wigner function of the system is calculated from the
density matrix $\rho_{\tau}$. We note that for experiments, there are various
well-developed methods for measuring Wigner function. For example, it can be
reconstructed from the homodyne measurements of the rotated quadrature phase
Vogel and Risken (1989) or direct measurements based on photon counting
Banaszek  _et al._ (1999). In the present scheme, this may be equivalently
performed on the emitted light from the microcavity. The Wigner function is
then sampled discretely on a grid in phase space and the samples arranged into
a vector, $\mathbf{f}$. A final signal is calculated as
$\bm{Y}^{\text{final}}=\bm{W}^{\text{out}}\bm{X}$, which is capable of input
classification. Here $\bm{W}^{\text{out}}$ is a trained output weight matrix
and the vector $\bm{X}=[1,\textbf{u},\textbf{f}]$, which is allowed to also
contain the input data pixels $u_{i}$ ($i=1,2,\cdots,16$) obtained as follows.
The input image with $28\times 28$ resolution is compressed into $4\times 4$
pixels (an equivalent procedure has been used in Ref. Ballarini  _et al._
(2020)), which are then converted into time-dependent signals with alternating
signs, such that each pixel corresponds to a different input signal value
$u_{i}$ in time, as shown in Fig. 1b. Note that even though the reduced
$4\times 4$ images are rather unrecognizable to humans, they can still be
recognized with the machine trained to do so. We show below that our method
can classify the $4\times 4$ resolution images and perform better than
classical polariton computers. Note that all the input signal values u are
weighted with a random mask (a $16\times 16$ matrix $\bm{W}^{\text{in}}$),
where its elements are randomly sampled from $0$ to
${W}^{\text{in}}_{\text{max}}$, before this image information enters the
quantum reservoir as time-dependent pumping amplitudes ($F$ and $P$) or phase
($\Theta$).

##  III III. Performance

Note that we have used $\Delta/\gamma=0.1$, $\alpha/\gamma=0.1$ and
$\gamma\tau/\hbar=1.2$ in our simulations. The $W^{\text{in}}_{\text{max}}$ is
optimised for pumping amplitudes $F/\gamma$, $P/\gamma$ and the phase
$\Theta$. The training procedure is performed on the output weights
$\bm{W}^{\text{out}}$ with ridge regression using a set of randomly chosen and
known images from the MNIST database. After the training, we keep the quantum
reservoir fixed and use a set of unknown images to test and evaluate its
efficiency for the image classification task. Details of the supervised
training of the polariton reservoir computing can be found in Appendix A.

We fix the number of training samples at $4000$ and testing samples at $1000$.
A separate analysis using a linear classifier method Lecun  _et al._ (1998),
results in a classification error rate of $29.8\%$ and $26.4\%$ for normalized
and unnormalized input signals u, respectively. Due to the compression of the
resolution of each image, it is difficult to achieve a low error rate with a
linear classifier. For our quantum polariton computing, however, considerably
lower error rate is possible. Our results for four different input encoding
methods ($\mathcal{SP},\mathcal{TP},\mathcal{SPTP},$ and
$\mathcal{SPTP}^{\prime}$) are shown in Fig. 2. Details of the influence of
physical parameters in the system are presented in Appendix B. To assess the
scaling of the performance, rather than considering scaling with the system
size, we consider scaling with a measure of the number of Fock states
available to the quantum reservoir. To do this, we introduce an effective
cutoff in the Fock space size, $N$. This effective cutoff $N$ naturally
appears for a finite pumping strength, which can only inject a finite number
of polaritons in the mode. Of course, in a physical system the Fock space is
unbounded (which can be accessed by increasing the pumping strength), but we
wish to emphasize that the scaling of Fock space size (which grows
quadratically with $N$) can be enhanced by the appropriate choice of
information encoding. For different pumping methods, the performance is indeed
found to improve from $\mathcal{SP}$ to $\mathcal{TP}$ to $\mathcal{SPTP}$,
and finally to $\mathcal{SPTP}^{\prime}$. This improvement originates from the
different ways the Wigner function contains the information (distribution and
squeezing) in regards to the different pumping scenarios, see details of the
analysis in Appendix C.

Figure 2: Classification error rate for different quantum reservoir sizes
(Fock space dimension) $N$ and different pumping methods: labelled
$\mathcal{SP}$ (single-photon); $\mathcal{TP}$ (two-photon); $\mathcal{SPTP}$
(single- and two-photon); and $\mathcal{SPTP}^{\prime}$ (single- and phase-
dependent two-photon) as written in the Hamiltonian of Eq. (2). Each data
point is obtained from the results of $10$ random realizations of the input
weights $\bm{W}^{\text{in}}$. Panels a and b present the lowest and average
error rate among the 10 realizations of the input weights, respectively, of
all the random realizations.

##  IV IV. Classical polariton reservoir

To demonstrate the impressive performance of the quantum polariton reservoir
in image recognition, we compare our system with classical polariton reservoir
computing. In this case, we consider a square polariton lattice with $N_{c}$
polariton modes having random nearest-neighbour coupling in a tight-binding
model. The dynamics is described by the general discrete nonlinear Schrödinger
equation:

$\displaystyle i\hbar\frac{\partial\psi_{n}}{\partial t}$ $\displaystyle=$
$\displaystyle
i\sum_{j=1,2}\bm{W}_{nj}^{\text{in}}\mathbf{v}_{j}+\sum\limits_{m\in\tilde{n}}\bm{K}_{nm}\psi_{m}+i\gamma\psi_{n}$
(4) $\displaystyle+\alpha|\psi_{n}|^{2}\psi_{n},$

where $\bm{W}^{\text{in}}$ (now a $N_{c}\times 2$ matrix) is an input mask
having random weights
$\bm{W}^{\text{in}}_{nj}\in[0,{W}^{\text{in}}_{\text{max}}]$, the vector
$\mathbf{v}=[1,u]$ containing the input signal, $\bm{K}_{nm}$ the random
nearest-neighbour coupling in the lattice, $\gamma$ the overall gain-loss
profile, $\alpha$ the Kerr nonlinear interaction strength and $\tilde{n}$
denotes the nearest neighbors of mode $n$. Note that in the classical
polariton reservoir computing, the obtained output from the reservoir is the
intensity of each mode, i.e., $|\psi_{n}|^{2}$, whereas in the quantum
polariton computing, we use the Wigner function. Similarly to the quantum
case, here $\bm{X}=[1,\textbf{u},\textbf{f}]$, where now the vector
$\textbf{f}=[|\psi_{1}|^{2},|\psi_{2}|^{2},\cdots,|\psi_{N_{c}}|^{2}]$.

In simulations, we keep the physical parameters the same as used in the
quantum scenario, i.e., $\alpha/\gamma=0.1$. As illustrated in Fig. 3, by
increasing the lattice size of the classical polariton reservoir, the error
rate saturates at $\sim 14\%$. To achieve the same (best) performance as in
the quantum polariton reservoir computing, one requires a classical reservoir
size of $100$ polariton modes for $\mathcal{SP}$ and $\mathcal{TP}$, and more
than $250$ polariton modes for $\mathcal{SPTP}$. For the lowest error rate
($13\%$) obtained from $\mathcal{SPTP}^{\prime}$, the classical polariton
computer still cannot be on par even with nearly $700$ polariton modes. We
note that the recent experiments in Ref. Ballarini  _et al._ (2020) reported
$86\%$ success rate for $4\times 4$ resolution with a 400-node classical
polariton reservoir. In our case, this is achieved theoretically with a single
quantum polariton mode.

Figure 3: Classification error rate (best and average performance) of the
classical polariton reservoir for different reservoir sizes $N_{c}$. Each data
point here is obtained from $10$ random realizations of the input weight
matrix $\bm{W}^{\text{in}}$. The four dashed lines correspond to the best
performance of the quantum polariton reservoir in Fig. 2.

##  V V. Quantum enhancement

By comparison with classical polariton reservoir computing, the quantum
polariton reservoir has a clear improvement in terms of performance and system
size. This is possible due to the quantum enhancement originating from the
high dimensionality of the polariton Fock space. As a measure of quantum
enhancement, the physical size of a quantum system is traditionally compared
with required size of a classical system to achieve the same level of
performance. In this respect, a classical reservoir computer requires hundreds
of polariton modes to achieve the same level of success rates shown by only a
single quantum polariton mode. This is clearly a superpolynomial quantum
enhancement. However, noting that a quantum polariton mode involves many Fock
states, another comparison would be between the physical size of a classical
reservoir and the effective number of qubits needed if one were to construct
the Fock space of the quantum mode. This comparison is typically used to
demonstrate quantum advantage in terms of information storage, see, e.g.,
Refs. Tacchino  _et al._ (2019); Marković and Grollier (2020).

The advantage in the quantum regime is illustrated in Fig. 4. The horizontal
axes in Fig. 4 are defined as the effective number of qubits
$N_{q}=\log_{2}N$, where $N$ is the effective Fock space dimension (minimum
number of Fock states required to describe the system). Ideally, the Fock
space dimension of a polariton mode is infinite. However, with a finite amount
of pumping strength only a finite number of Fock states will be occupied,
which defines a finite effective dimension. Each data point in Fig. 4
corresponds to a specific error rate for image recognition; for example, the
point indicated by the arrow in panel $\mathbf{b}$ represents achieving an
error rate $\leq$ $14\%$, in which case it requires less than four effective
qubits for the quantum polariton reservoir and $500$ modes for the classical
polariton scenario. It is apparent from Fig. 4 that the quantum polariton
reservoir showcases superpolynomial enhanced performance compared to the
classical polariton reservoir. Note that the number of classical modes $N_{c}$
is plotted in logarithmic scale in both panels.

Figure 4: Superpolynomial enhancement of the quantum polariton reservoir
computing. The vertical axis in each panel represents the number of polariton
modes required for achieving a specific error rate, while the horizontal axis
is the number of effective qubits required in the quantum reservoir computing.
Panels a and b represent the best and average performance, respectively, of
the quantum and classical reservoir.

##  VI VI. Conclusion

In this paper, a quantum reservoir network that is formed by the Fock space of
quantum polaritons is shown to possess the capability for effective image
classification task. By utilizing different pumping methods, i.e., single-
photon and phase-dependent two-photon pumping, we show performances with an
error rate as low as $13\%$ with practical parameters (e.g.,
$\alpha/\gamma=0.1$). Our results show that by exploiting the Fock space
dimension of one quantum polariton mode, one can achieve lower error rate for
image classification task, than that obtained from a classical reservoir with
$\sim 700$ polariton modes. Finally, by comparing the required system sizes of
the quantum polariton reservoir (effective number of qubits) and the classical
polariton case for achieving the same classification efficiency,
superpolynomial quantum enhancement is demonstrated. Our results motivate
applications exploiting the quantum nature of exciton-polaritons.

##  Appendix A APPENDIX A: Supervised Training

In the supervised training, we offer the reservoir network a set of known
images from the MNIST database, which correspond to the input signals
$\mathbf{u}_{k}$ with the subscript denoting the $k$th image, where each
$\mathbf{u}_{k}=[u_{1},u_{2},\cdots,u_{16}]$ and $k=1,2,\cdots,N_{tr}$. Also,
before the signals enter the reservoir, we apply a random mask
$\bm{W}^{\text{in}}$ to each of the input signals $\mathbf{u}_{k}$. With the
known images, we have the corresponding desired target outputs
$\mathbf{Y}^{\text{target}}_{k}$. The output we obtain directly from the
reservoir network is its Wigner function, arranged in a vector form
$\mathbf{f}_{k}$. We take the output layer of the reservoir as a vector
$\mathbf{X}_{k}=[1,\mathbf{u}_{k},\mathbf{f}_{k}]$ and the relation between
the final output and the output layer is given by
$\mathbf{Y}^{\text{final}}_{k}=\bm{W}^{\text{out}}\mathbf{X}_{k}$. For a set
of known examples, we have

$\displaystyle\mathcal{Y}^{\text{final}}=\bm{W}^{\text{out}}\mathcal{X},$ (A1)

where
$\mathcal{Y}^{\text{final}}=[\mathbf{Y}_{1}^{\text{final}};\mathbf{Y}_{2}^{\text{final}};...;\mathbf{Y}_{N_{tr}}^{\text{final}}]$
is the collection of all the final outputs and
$\mathcal{X}=[\mathbf{X}_{1};\mathbf{X}_{2};...;\mathbf{X}_{N_{tr}}]$ is the
corresponding collection of the obtained outputs. The purpose of the training
is to minimize the difference between $\mathcal{Y}^{\text{final}}$ and
$\mathcal{Y}^{\text{target}}$, which stands for the collection of all the
target outputs, i.e.,
$\mathcal{Y}^{\text{target}}=[\mathbf{Y}_{1}^{\text{target}};\mathbf{Y}_{2}^{\text{target}};...;\mathbf{Y}_{N_{tr}}^{\text{target}}]$.
We apply the ridge regression Hoerl and Kennard (1970) here to get the
optimized output weight matrix $\bm{W}^{\text{out}}$:

$\displaystyle\bm{W}^{\text{out}}=\mathcal{Y}^{\text{target}}{\mathcal{X}}^{T}(\mathcal{X}\mathcal{X}^{T}+{\beta}\bm{I})^{-1},$
(A2)

where $\mathcal{X}^{T}$ is the transpose of the matrix $\mathcal{X}$, $I$ is
an identity matrix and $\beta$ is the regularization coefficient. The term
$\beta I$ is used to control the $\bm{W}^{\text{out}}$, in case there is small
fluctuation in the input signals, the results still converge.

##  Appendix B APPENDIX B: Analysis of Physical Parameters

Before the input signals ($\mathbf{u}$) enter the reservoir network, a random
mask $\bm{W}^{\text{in}}$ is applied to it. The elements in the mask are
randomly distributed from $0$ to ${W}_{\text{max}}^{\text{in}}$. This way, the
actual input signals entering the reservoir will be
$\bm{W}^{\text{in}}\mathbf{u}$. In other words, the amplitude of the input
mask determines the pumping strength and/or phase amplitude. We take the
$\mathcal{SP}$ pumping method as an example, vary the
${W}^{\text{in}}_{\text{max}}$ and analyse the classification error rate. As
shown in Fig. B1, with increasing of ${W}_{\text{max}}^{\text{in}}$, the
classification error rate goes down.

Figure B1: Average and lowest classification error rate for different input
weight amplitudes (${W}^{\text{in}}_{\text{max}}$), for the $\mathcal{SP}$
pumping scheme. Each point is obtained from $10$ random realizations of the
matrix $\bm{W}^{\text{in}}$. $N$ is fixed at 10.

The significance of nonlinearity in the system is analysed as well for the
$\mathcal{SP}$ pumping method. To make the role of nonlinearity more
prominent, we vary $\alpha$ while fixing the randomly chosen mask
$\bm{W}^{\text{in}}$ and $\gamma$. As can be seen in Fig. B2, with increasing
the nonlinearity, the classification error rate drops, which indicates the
quantum polariton reservoir works better in stronger nonlinear regime.

Figure B2: Effect of nonlinearity in the system, showing the classification
error rate under different nonlinear strength.

In real situations, there might be a fluctuation in the detuning $\Delta$
coming from the fluctuation in the driving laser frequency $\omega_{p}$. To
investigate this, we added a random time-dependent energy fluctuation term
$\eta(t)\in[-1,1]\delta$ (where $\delta$ is the fluctuation strength) in the
detuning. As an example, we plotted the error rate against the fluctuation
strength for the single photon pumping case ($\mathcal{SP}$) in Fig. B3. It
can be seen that the classification error rate is almost unchanged ($\sim
24.5\%$), which holds true even if the fluctuation-to-linewidth ratio is
$\delta/\gamma\sim 0.1$, corresponding to $\delta=0.05$ meV for
$\gamma=0.5$meV..

We would like to note that a frequency stabilized laser may be used as a
reference for highly stable pump. In this case, for example, a fractional
frequency stability $\sim 5\times 10^{-16}/\sqrt{\tau_{m}}$ ($\tau_{m}$ is
measurement time in seconds) has been demonstrated in Ref. Jiang  _et al._
(2011). For polaritons, a typical laser energy is $\sim 2$ eV, which gives a
fractional frequency variation $2.5\times 10^{-5}$ in our case. In the case of
highly stable lasers, for $\tau_{m}=0.1$ ps, the fractional frequency
stability is $\sim 1.6\times 10^{-9}$, which is much smaller than what our
system requires.

Figure B3: Dependence of classification error rate on the fluctuation strength
of the driving frequency for the single photon pumping case.

The two photon pumping can be realized with an inverse four-wave mixing
process, e.g., see Ref. Kyriienko and Liew (2014). This theoretically proposed
configuration is similar to routinely studied four-wave mixing in
microcavities Saba  _et al._ (2001); Wu  _et al._ (2021); Baumberg  _et al._
(2000); Sanvitto  _et al._ (2006); Savvidis  _et al._ (2000); Kundermann  _et
al._ (2003), only instead of considering a polariton amplifier with
simultaneous laser excitation of signal and pump modes, one considers
simultaneous excitation of signal and idler modes, which results in effective
two particle excitation of the pump mode. In such case
$P=\beta\sqrt{n_{s}n_{i}}$, where $\beta$ is the polariton-polariton
interaction strength, $n_{s}$ is the density of polaritons in the signal
state, and $n_{i}$ is the density in the idler state. As $\sqrt{n_{s}n_{i}}$
is effectively a measure of the density in both signal and idler states (i.e.,
the geometric mean), a typical value of $P$ would correspond to the typical
blueshifts measured in experiments under resonant excitation. Indeed these can
be anywhere between 0 and 1 meV. Varying this value is equivalent to changing
the strength of laser fields, which is routinely done, potentially with an
error as small as $0.02\%$, see for example Ref. Abbaspour  _et al._ (2015).
In our case, the smallest change of P is $\sim 0.2$ meV, which is much higher
than $2\times 10^{-4}$ meV mentioned above.

Last, we also note that the total time should be within the polariton
lifetime. Our simulations indicate that if we go beyond the polariton
lifetime, the error rate increases. This is intuitive as if we go beyond the
polariton lifetime, the information from the earlier input signals fed into
the system is already lost before the later signals are implemented. This is a
key difference from a classical system, where an incoherent pump could
compensate the loss due to finite lifetime. In a quantum reservoir, this is
not possible as incoherently injected polaritons will not have the same
quantum correlation with rest of the polaritons injected by the input signal.

Figure C1: Distribution of Wigner function under different pumping methods.
Panel a-d correspond to the Wigner function with $\mathcal{SP}$,
$\mathcal{TP}$, $\mathcal{SPTP}$ and $\mathcal{SPTP}^{\prime}$ pumping
methods, respectively. Panels e-h are the corresponding ones with higher
pumping strength.

##  Appendix C APPENDIX C: Wigner Function of Different Pumping Methods

As shown in the main text, we have a clear improvement by utilising different
pumping methods for the quantum reservoir network, the origin of which we now
analyse. As we use Wigner function as the output from the quantum reservoir
network, we shall check its distribution for each of the pumping schemes and
strength. Fig. C1a, b, c and d show the quasiprobability distribution in phase
space under weaker pumping strength for the scenario of $\mathcal{SP}$,
$\mathcal{TP}$, $\mathcal{SPTP}$, and $\mathcal{SPTP}^{\prime}$, respectively
(see Eq. (2) in the main text and the explanation that follows). Furthermore,
by increasing the pumping strength, the corresponding Wigner functions appear
as in Fig. C1e-h.

It is clear from the $\mathcal{SP}$ pumping scheme in Fig. C1 panels a and e,
that by increasing the pumping strength, the Wigner function retains the
Gaussian-shape, but the center moves. In this case, for different input
signals, the center of quasiprobability distribution in phase space is driven
to different positions. A different mechanism happens for the $\mathcal{TP}$
pumping method, see Fig. C1b and f. For different pumping strength, the
distribution stays in the center but they have different “squeezed” shapes.
Hence, for $\mathcal{SP}$ and $\mathcal{TP}$ one can think that, in the phase
space, we have one degree of freedom for each method, which is position
($\mathcal{SP}$) and squeezing ($\mathcal{TP}$). To some extent, this explains
why the error rates of $\mathcal{SP}$ and $\mathcal{TP}$ pumping methods are
quite close.

However, if we look at Fig. C1c and g, with stronger pumping strength, the
position and the squeezing of the distribution both change, which means that
under the combined pumping schemes $\mathcal{SPTP}$, we have two degrees of
freedom. With this scheme, different input signals will provide us with more
varied outputs, and this improves the classification process. We go further to
the fourth pumping method ($\mathcal{SPTP}^{\prime}$), in which we also encode
the input information with a phase $\Theta$. As shown in Fig. C1d and h, when
we increase the pumping strength, we obtained a Wigner function with different
center position, squeezing strength and direction. With these three degrees of
freedom, we are able to get a classification error rate as low as $13\%$.

The Wigner functions shown in Fig. C1 are not negative (which is normally used
to indicate the presence of quantumness). We note that quantum enhancement in
this paper is associated with the exploitation of the largely available Fock
space as compared to the classical multi-network polariton computer. The
different pumping schemes then serve as different “encoding” methods in which
information is fed into the largely available space in the quantum case.

## Acknowledgement

This work was supported by the Singapore Ministry of Education under its AcRF
Tier 2 grant MOE2019-T2-1-004.

## References

  * Carusotto and Ciuti (2013) I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
  * Kasprzak  _et al._ (2006) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Nature 443, 409 (2006).
  * Sanvitto and Kéna-Cohen (2016) D. Sanvitto and S. Kéna-Cohen, Nat. Mater. 15, 1061 (2016).
  * Fraser (2017) M. D. Fraser, Semicond. Sci. Technol. 32, 093003 (2017).
  * Leyder  _et al._ (2007) C. Leyder, T. C. H. Liew, A. V. Kavokin, I. A. Shelykh, M. Romanelli, J. P. Karr, E. Giacobino, and A. Bramati, Phys. Rev. Lett. 99, 196402 (2007).
  * Cancellieri  _et al._ (2015) E. Cancellieri, J. K. Chana, M. Sich, D. N. Krizhanovskii, M. S. Skolnick, and D. M. Whittaker, Phys. Rev. B 92, 174528 (2015).
  * Nguyen  _et al._ (2013) H. S. Nguyen, D. Vishnevsky, C. Sturm, D. Tanese, D. Solnyshkov, E. Galopin, A. Lemaître, I. Sagnes, A. Amo, G. Malpuech, and J. Bloch, Phys. Rev. Lett. 110, 236601 (2013).
  * Niemietz  _et al._ (2016) D. Niemietz, J. Schmutzler, P. Lewandowski, K. Winkler, M. Aßmann, S. Schumacher, S. Brodbeck, M. Kamp, C. Schneider, S. Höfling, and M. Bayer, Phys. Rev. B 93, 235301 (2016).
  * Gao  _et al._ (2012) T. Gao, P. S. Eldridge, T. C. H. Liew, S. I. Tsintzos, G. Stavrinidis, G. Deligeorgis, Z. Hatzopoulos, and P. G. Savvidis, Phys. Rev. B 85, 235102 (2012).
  * De Giorgi  _et al._ (2012) M. De Giorgi, D. Ballarini, E. Cancellieri, F. M. Marchetti, M. H. Szymanska, C. Tejedor, R. Cingolani, E. Giacobino, A. Bramati, G. Gigli, and D. Sanvitto, Phys. Rev. Lett. 109, 266407 (2012).
  * Dreismann  _et al._ (2016) A. Dreismann, H. Ohadi, Y. del Valle-Inclan Redondo, R. Balili, Y. G. Rubo, S. I. Tsintzos, G. Deligeorgis, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Nat. Mater. 15, 1074 (2016).
  * Grosso  _et al._ (2014) G. Grosso, S. Trebaol, M. Wouters, F. Morier-Genoud, M. T. Portella-Oberli, and B. Deveaud, Phys. Rev. B 90, 045307 (2014).
  * Rozas  _et al._ (2020) E. Rozas, J. Beierlein, A. Yulin, M. Klaas, H. Suchomel, O. Egorov, I. A. Shelykh, U. Peschel, C. Schneider, S. Klembt, S. Höfling, M. D. Martín, and L. Viña, Adv. Opt. Mater. 8, 2070072 (2020).
  * Klaas  _et al._ (2019) M. Klaas, J. Beierlein, E. Rozas, S. Klembt, H. Suchomel, T. H. Harder, K. Winkler, M. Emmerling, H. Flayac, M. D. Martín, L. Viña, S. Höfling, and C. Schneider, Appl. Phys. Lett. 114, 061102 (2019).
  * Winkler  _et al._ (2017) K. Winkler, H. Flayac, S. Klembt, A. Schade, D. Nevinskiy, M. Kamp, C. Schneider, and S. Höfling, Phys. Rev. B 95, 201302 (2017).
  * Ballarini  _et al._ (2013) D. Ballarini, M. De Giorgi, E. Cancellieri, R. Houdré, E. Giacobino, R. Cingolani, A. Bramati, G. Gigli, and D. Sanvitto, Nat. Commun. 4, 1778 (2013).
  * Lewandowski  _et al._ (2017) P. Lewandowski, S. M. H. Luk, C. K. P. Chan, P. T. Leung, N. H. Kwong, R. Binder, and S. Schumacher, Opt. Express 25, 31056 (2017).
  * Zasedatelev  _et al._ (2019) A. V. Zasedatelev, A. V. Baranikov, D. Urbonas, F. Scafirimuto, U. Scherf, T. Stöferle, R. F. Mahrt, and P. G. Lagoudakis, Nat. Photonics 13, 378 (2019).
  * Flayac and Savenko (2013) H. Flayac and I. G. Savenko, Appl. Phys. Lett. 103, 201105 (2013).
  * Marsault  _et al._ (2015) F. Marsault, H. S. Nguyen, D. Tanese, A. Lemaître, E. Galopin, I. Sagnes, A. Amo, and J. Bloch, Appl. Phys. Lett. 107, 201115 (2015).
  * Lerario  _et al._ (2017) G. Lerario, D. Ballarini, A. Fieramosca, A. Cannavale, A. Genco, F. Mangione, S. Gambino, L. Dominici, M. De Giorgi, G. Gigli, and D. Sanvitto, Light Sci. Appl. 6, e16212 (2017).
  * Cerna  _et al._ (2013) R. Cerna, Y. Léger, T. K. Paraïso, M. Wouters, F. Morier-Genoud, M. T. Portella-Oberli, and B. Deveaud, Nat. Commun. 4, 2008 (2013).
  * Kalinin and Berloff (2019) K. P. Kalinin and N. G. Berloff, Phys. Rev. B 100, 245306 (2019).
  * Kalinin and Berloff (2018) K. P. Kalinin and N. G. Berloff, Sci. Rep. 8, 17791 (2018).
  * Banerjee and Liew (2020) R. Banerjee and T. C. H. Liew, New J. Phys. 22, 103062 (2020).
  * Espinosa-Ortega and Liew (2015) T. Espinosa-Ortega and T. C. H. Liew, Phys. Rev. Lett. 114, 118101 (2015).
  * Ohadi  _et al._ (2017) H. Ohadi, A. J. Ramsay, H. Sigurdsson, Y. del Valle-Inclan Redondo, S. I. Tsintzos, Z. Hatzopoulos, T. C. H. Liew, I. A. Shelykh, Y. G. Rubo, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. Lett. 119, 067401 (2017).
  * Töpfer  _et al._ (2020) J. D. Töpfer, I. Chatzopoulos, H. Sigurdsson, T. Cookson, Y. G. Rubo, and P. G. Lagoudakis, arXiv e-prints , arXiv:2007.06690 (2020).
  * Tanaka  _et al._ (2019) G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, Neural networks : the official journal of the International Neural Network Society 115, 100 (2019).
  * Jaeger and Haas (2004) H. Jaeger and H. Haas, Science 304, 78 (2004).
  * Lukoševičius and Jaeger (2009) M. Lukoševičius and H. Jaeger, Comput. Sci. Rev. 3, 127 (2009).
  * Enel  _et al._ (2016) P. Enel, E. Procyk, R. Quilodran, and P. F. Dominey, PLoS Comput. Biol. 12, e1004967 (2016).
  * Lukoševičius (2012) M. Lukoševičius, “A practical guide to applying echo state networks,” in  _Neural Networks: Tricks of the Trade: Second Edition_, edited by G. Montavon, G. B. Orr, and K.-R. Müller (Springer Berlin Heidelberg, 2012) p. 659.
  * Opala  _et al._ (2019) A. Opala, S. Ghosh, T. C. H. Liew, and M. Matuszewski, Phys. Rev. Appl. 11, 064029 (2019).
  * Xu  _et al._ (2020) H. Xu, S. Ghosh, M. Matuszewski, and T. C. H. Liew, Phys. Rev. Appl. 13, 064074 (2020).
  * Ballarini  _et al._ (2020) D. Ballarini, A. Gianfrate, R. Panico, A. Opala, S. Ghosh, L. Dominici, V. Ardizzone, M. De Giorgi, G. Lerario, G. Gigli, T. C. H. Liew, M. Matuszewski, and D. Sanvitto, Nano Lett. 20, 3506 (2020).
  * Mirek  _et al._ (2021) R. Mirek, A. Opala, P. Comaron, M. Furman, M. Król, K. Tyszka, B. Seredyński, D. Ballarini, D. Sanvitto, T. C. H. Liew, W. Pacuski, J. Suffczyński, J. Szczytko, M. Matuszewski, and B. Piętka, Nano Letters (2021), 10.1021/acs.nanolett.0c04696.
  * Savasta  _et al._ (2005) S. Savasta, O. D. Stefano, V. Savona, and W. Langbein, Phys. Rev. Lett. 94, 246401 (2005).
  * Karr  _et al._ (2004) J. P. Karr, A. Baas, R. Houdré, and E. Giacobino, Phys. Rev. A 69, 031802 (2004).
  * Boulier  _et al._ (2014) T. Boulier, M. Bamba, A. Amo, C. Adrados, A. Lemaitre, E. Galopin, I. Sagnes, J. Bloch, C. Ciuti, E. Giacobino, and A. Bramati, Nat. Commun. 5, 3260 (2014).
  * Cuevas  _et al._ (2018) Á. Cuevas, J. C. López Carreño, B. Silva, M. De Giorgi, D. G. Suárez-Forero, C. Sánchez Muñoz, A. Fieramosca, F. Cardano, L. Marrucci, V. Tasco, G. Biasiol, E. del Valle, L. Dominici, D. Ballarini, G. Gigli, P. Mataloni, F. P. Laussy, F. Sciarrino, and D. Sanvitto, Sci. Adv. 4, eaao6814 (2018).
  * Suárez-Forero  _et al._ (2020) D. G. Suárez-Forero, V. Ardizzone, S. F. C. da Silva, M. Reindl, A. Fieramosca, L. Polimeno, M. De Giorgi, L. Dominici, L. N. Pfeiffer, G. Gigli, _et al._ , Light Sci. Appl. 9, 1 (2020).
  * Muñoz-Matutano  _et al._ (2019) G. Muñoz-Matutano, A. Wood, M. Johnsson, X. Vidal, B. Q. Baragiola, A. Reinhard, A. Lemaître, J. Bloch, A. Amo, G. Nogues, _et al._ , Nat. Mater. 18, 213 (2019).
  * Delteil  _et al._ (2019) A. Delteil, T. Fink, A. Schade, S. Höfling, C. Schneider, and A. İmamoğlu, Nat. Mater. 18, 219 (2019).
  * Verger  _et al._ (2006) A. Verger, C. Ciuti, and I. Carusotto, Phys. Rev. B 73, 193306 (2006).
  * Kyriienko  _et al._ (2020) O. Kyriienko, D. N. Krizhanovskii, and I. A. Shelykh, Phys. Rev. Lett. 125, 197402 (2020).
  * Rosenberg  _et al._ (2018) I. Rosenberg, D. Liran, Y. Mazuz-Harpaz, K. West, L. Pfeiffer, and R. Rapaport, Sci. Adv. 4, eaat8880 (2018).
  * Jia  _et al._ (2018) N. Jia, N. Schine, A. Georgakopoulos, A. Ryou, L. W. Clark, A. Sommer, and J. Simon, Nat. Phys. 14, 550 (2018).
  * Stefanatos and Paspalakis (2020) D. Stefanatos and E. Paspalakis, Phys. Rev. A 102, 013716 (2020).
  * Fujii and Nakajima (2017) K. Fujii and K. Nakajima, Phys. Rev. Appl. 8, 024030 (2017).
  * Nakajima  _et al._ (2019) K. Nakajima, K. Fujii, M. Negoro, K. Mitarai, and M. Kitagawa, Phys. Rev. Applied 11, 034021 (2019).
  * Chen and Nurdin (2019) J. Chen and H. I. Nurdin, Quantum Inf. Process. 18, 198 (2019).
  * Chen  _et al._ (2020) J. Chen, H. I. Nurdin, and N. Yamamoto, Phys. Rev. Applied 14, 024065 (2020).
  * Govia  _et al._ (2021) L. C. G. Govia, G. J. Ribeill, G. E. Rowlands, H. K. Krovi, and T. A. Ohki, Physical Review Research 3, 013077 (2021).
  * Kalfus  _et al._ (2021) W. D. Kalfus, G. J. Ribeill, G. E. Rowlands, H. K. Krovi, T. A. Ohki, and L. C. G. Govia, arXiv e-prints , arXiv:2101.11729 (2021).
  * Olivares (2012) S. Olivares, Eur. Phys. J. Spec. Top. 203, 3 (2012).
  * Bartolo  _et al._ (2016) N. Bartolo, F. Minganti, W. Casteels, and C. Ciuti, Phys. Rev. A 94, 033841 (2016).
  * Vogel and Risken (1989) K. Vogel and H. Risken, Phys. Rev. A 40, 2847 (1989).
  * Banaszek  _et al._ (1999) K. Banaszek, C. Radzewicz, K. Wódkiewicz, and J. S. Krasiński, Phys. Rev. A 60, 674 (1999).
  * Lecun  _et al._ (1998) Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner, Proc. IEEE 86, 2278 (1998).
  * Tacchino  _et al._ (2019) F. Tacchino, C. Macchiavello, D. Gerace, and D. Bajoni, npj Quantum Inf. 5, 26 (2019).
  * Marković and Grollier (2020) D. Marković and J. Grollier, Appl. Phys. Lett. 117, 150501 (2020).
  * Hoerl and Kennard (1970) A. E. Hoerl and R. W. Kennard, Technometrics 12, 55 (1970).
  * Jiang  _et al._ (2011) Y. Y. Jiang, A. D. Ludlow, N. D. Lemke, R. W. Fox, J. A. Sherman, L.-S. Ma, and C. W. Oates, Nat. Photonics 5, 158 (2011).
  * Kyriienko and Liew (2014) O. Kyriienko and T. C. H. Liew, Phys. Rev. A 90, 063805 (2014).
  * Saba  _et al._ (2001) M. Saba, C. Ciuti, J. Bloch, V. Thierry-Mieg, R. André, L. S. Dang, S. Kundermann, A. Mura, G. Bongiovanni, J. L. Staehli, and B. Deveaud, Nature 414, 731 (2001).
  * Wu  _et al._ (2021) J. Wu, S. Ghosh, R. Su, A. Fieramosca, T. C. H. Liew, and Q. Xiong, Nano Letters 21, 3120 (2021).
  * Baumberg  _et al._ (2000) J. J. Baumberg, P. G. Savvidis, R. M. Stevenson, A. I. Tartakovskii, M. S. Skolnick, D. M. Whittaker, and J. S. Roberts, Phys. Rev. B 62, R16247 (2000).
  * Sanvitto  _et al._ (2006) D. Sanvitto, D. N. Krizhanovskii, D. M. Whittaker, S. Ceccarelli, M. S. Skolnick, and J. S. Roberts, Phys. Rev. B 73, 241308 (2006).
  * Savvidis  _et al._ (2000) P. G. Savvidis, J. J. Baumberg, R. M. Stevenson, M. S. Skolnick, D. M. Whittaker, and J. S. Roberts, Phys. Rev. Lett. 84, 1547 (2000).
  * Kundermann  _et al._ (2003) S. Kundermann, M. Saba, C. Ciuti, T. Guillet, U. Oesterle, J. L. Staehli, and B. Deveaud, Phys. Rev. Lett. 91, 107402 (2003).
  * Abbaspour  _et al._ (2015) H. Abbaspour, G. Sallen, S. Trebaol, F. Morier-Genoud, M. T. Portella-Oberli, and B. Deveaud, Phys. Rev. B 92, 165303 (2015).


# Finding the Optimum Design of Large Gas Engines Prechambers Using CFD and
Bayesian Optimization

Stefan Posch  LEC GmbH Clemens Gößnitzer  LEC GmbH Franz Rohrhofer  Know-
Center GmbH  Bernhard C. Geiger  Know-Center GmbH  Andreas Wimmer  LEC GmbH
Institut für Thermodynamik und nachhaltige Antriebssysteme, TU Graz

###### ?abstractname?

The turbulent jet ignition concept using prechambers is a promising solution
to achieve stable combustion at lean conditions in large gas engines, leading
to high efficiency at low emission levels. Due to the wide range of design and
operating parameters for large gas engine prechambers, the preferred method
for evaluating different designs is computational fluid dynamics (CFD), as
testing in test bed measurement campaigns is time-consuming and expensive.
However, the significant computational time required for detailed CFD
simulations due to the complexity of solving the underlying physics also
limits its applicability. In optimization settings similar to the present
case, i.e., where the evaluation of the objective function(s) is
computationally costly, Bayesian optimization has largely replaced classical
design-of-experiment. Thus, the present study deals with the computationally
efficient Bayesian optimization of large gas engine prechambers design using
CFD simulation. Reynolds-averaged-Navier-Stokes simulations are used to
determine the target values as a function of the selected prechamber design
parameters. The results indicate that the chosen strategy is effective to find
a prechamber design that achieves the desired target values.

##  1 Introduction

Internal combustion engines (ICE) will continue to play an important role in
the generation of energy and the transportation of goods. The high energy
density of the used fuels in combination with the robust and reliable
technology, the ICE will still be the first choice for applications in which
pure electrified systems make little sense as it is the case for the shipping
of goods. Furthermore, ICE will be a key part for future energy management by
the reconversion of chemical energy storage media such as hydrogen into
electricity or to ensure grid stability which is induced by the increasing use
of alternative energy sources such as wind and solar energy [GBV+19]. Hence,
the research on increasing engine efficiencies under the fulfilment of
stringent emission legislation will continue regardless of whether classical
carbon-based or alternative carbon-free fuels are used. In the case of spark
ignited (SI) engines, a common approach is the operation under lean conditions
which increases the engine efficiency by increasing the specific heat capacity
of the working fluid, reducing heat losses through the cylinder walls, and
reducing pumping losses in part load operation [ZALY22]. However, the
operation at lean conditions comes with the drawback of unstable combustion
and thus, high cyclic fluctuations that can even lead to misfiring cycles. The
enhancement of the ignition energy by increasing the performance of the
ignition system helps to overcome these issues. A well-established concept to
increase the ignition efficiency of lean mixture SI engines is the turbulent
jet ignition (TJI) by the use of a separate, small (compared to the main
combustion chamber) volume known as prechamber. The mixture is ignited by the
spark plug in the prechamber and the flame propagates through the volume
increasing the pressure and force hot products and active radicals to flow
into the main combustion chamber via overflow bores. Thus, the ignition of the
mixture in the main combustion chamber is induced at multiple locations due to
the penetrating turbulent jets coming from the prechamber. However, the design
of a prechamber and the choice of the appropriate operating conditions to
achieve a highly efficient TJI system in combination with the main combustion
chamber is not straightforward due to the large number of parameters and thus
degrees of freedom. Therefore, experimental investigations as shown by
Roethlisberger and Favrat [RF03] or Novella et al. [NPGS+20] are complex as
well as cost intense and thus limit the number of design variants that can be
tested. To decrease the experimental effort, the state-of-the-art procedure is
to apply numerical simulation, i.e. computational fluid dynamics (CFD), to
increase the number of design and operation parameter variants by virtual
testing and to achieve a comprehensive view on the physical effects which
measurements are not able to provide.

Various studies regarding the application of CFD to evaluate prechamber
efficiencies and related values can be found in literature. Benajes et al.
[BNGS+20] investigated the scavenging and turbulence distributions of several
prechamber designs under ultra-lean conditions. The results of the study
indicate that the reduction of the laminar flame speed due to the ultra-lean
conditions significantly decreases the quality of the turbulent jets exiting
the prechamber. Silva et al. [SSH+20] carried out CFD combustion simulations
to investigate the influence of several design parameters on the engine
combustion characteristics. The authors stated that the diameter of the
overflow bores impacts the peak pressure and the residence time of the main
charge in the prechamber. Furthermore, the results show a large influence of
the neck diameter on the pressure increase in the prechamber as well es the
aforementioned residual time. Related works can be found in Feng et al.
[FZQ+18], Winter et al. [WSP+20], or Zhang et al. [ZYE+20]. In addition to
common CFD simulation, data-driven approaches and machine learning (ML)
concepts are increasingly finding their way into the design optimization of
ICE components to facilitate the discovery of appropriate prechamber
configurations. Ge et al. [GBY+21] compared different ML models to develop a
surrogate model which represents combustion CFD simulations. The trained model
is further used in a genetic algorithm to find an optimum prechamber design.
To decrease the number of required CFD simulations, a Bayesian updating
optimization method was applied. In a similar study, Silva et al. [SMB+22]
compared design of experiments (DOE) optimization with different ML surrogate
models and combinations of optimization techniques via a so-called automated
super-learner framework. The results indicate that the proposed method
performance is superior compared to the classical DOE approach. Although these
studies use ML surrogate models for the optimization process, the training of
the models requires a certain number of combustion CFD simulations which are
computationally expensive. To overcome this issue, Posch et al. [PWZ+21]
presented a method to use condensed CFD simulations to train ML models, which
are able to predict trends in the prechamber performance, thus helping to
reduce the computational effort. Following the aforementioned studies of Ge et
al. [GBY+21] and Silva et al. [SMB+22], the present study deals with the
optimization of a prechamber design by the use of CFD and appropriate
optimization techniques. A Bayesian optimization (BO) algorithm based on
Gaussian processes is applied which is initialized by a DOE-based CFD data
set. To reduce the computational effort, the CFD simulation only covers the
range from bottom dead center until ignition timing (IT) rather than the
entire cycle. This simplification can be justified as follows: (i) The present
study focuses solely on the prechamber behaviour and thus gas exchange through
the valves in the prechamber is neglected. (ii) Since the flow field around
the spark plug at IT has a significant influence on the flame kernel
development [GKK+19, GG21, PGO+22] and furthermore on the combustion in the
pre- and main chamber, the optimization process focuses on the turbulent
kinetic energy and velocity magnitude around the spark plug at IT. The results
show that the proposed method is able to find a prechamber design that covers
three design parameters and that fulfills the target values in terms of
maximum turbulent kinetic energy by considering a certain maximum velocity
magnitude around the spark plug. Figure 1 shows the base prechamber design
with the three input parameters for the optimization process.

$d_{bottle}$

$h_{neck}$

$d_{bore}$

?figurename? 1: Base prechamber design and the three input parameters

##  2 Methodology

###  2.1 Engine specification

The engine considered here was adapted from an existing engine. The main
chamber was simplified to a perfect cylinder, and the prechamber had variable
geometrical sizes as described later. The simplified engine had a bore of
$190\text{\,}\mathrm{m}\mathrm{m}$, a stroke of
$220\text{\,}\mathrm{m}\mathrm{m}$, a displaced volume of about
$6.2\text{\,}\mathrm{L}$ and a speed of
$1500\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$. No valves or gas exchange were
modelled. Thus, the time of interest for this study covered the engine cycle
ranging from bottom dead centre at
$540\text{\,}\mathrm{C}\mathrm{A}\mathrm{D}$ until ignition timing at
$695\text{\,}\mathrm{C}\mathrm{A}\mathrm{D}$.

###  2.2 CFD setup

The Reynolds-averaged Navier-Stokes (RANS) CFD simulations were performed with
OpenFOAM (version from November 2022 [Opea] from The OpenFOAM Foundation
[Opeb]). Meshing was done with snappyHexMesh. The grid size was
$4\text{\,}\mathrm{m}\mathrm{m}$ in the cylinder,
$0.25\text{\,}\mathrm{m}\mathrm{m}$ in the prechamber and
$0.125\text{\,}\mathrm{m}\mathrm{m}$ in the prechamber bores and near the
spark plug, leading to approximately 1 mio cells. Turbulence was modelled by
the RANS $k$-$\omega$ SST two-equation model. Initial temperature and pressure
were set to $400\text{\,}\mathrm{K}$ and
$3\text{\,}\mathrm{b}\mathrm{a}\mathrm{r}$, respectively. The time step was
adapted to keep the Courant number below $2$. The PISO algorithm with three
correctors was used for pressure-velocity coupling. Additionally, two non-
orthogonal correctors were used to improve convergence of the pressure
prediction. Furthermore, second-order bounded discretization schemes were used
for gradient and divergence, whereas a first-order implicit time
discretization was used. Since the maximum observed speeds exceed a Mach
number of $0.3$, the transonic option was set for all simulations. Mesh
movement was done with the standard OpenFOAM layered engine approach. Mass
conservation was monitored by integration of the density over the volume, to
check if the total mass of the closed volume remains constant, which indeed
was the case for all simulations. A more detailed description of the
simulation setup can be found in [GPP22]. For the evaluation of the quantities
near the spark plug, all cells within a sphere with radius of
$0.5\text{\,}\mathrm{m}\mathrm{m}$ just outside of the spark plug gap were
monitored. The relevant simulation output, namely the turbulence intensity
measured by the volume-average value of $k$ and velocity magnitude $|\vec{v}|$
are both evaluated inside this sphere.

###  2.3 Optimization strategy

Three geometrical parameters of the prechamber were adjusted in this study to
find the optimal shape: the diameter of the prechamber bottle
($d_{\mathrm{bottle}}$), the diameter of the overflow bores
($d_{\mathrm{bore}}$), and the height of the prechamber neck
($h_{\mathrm{neck}}$). The optimization ranges were chosen as
$d_{\mathrm{bottle}}\in[8,12]$, $d_{\mathrm{bore}}\in[0.75,1.15]$ and
$h_{\mathrm{neck}}\in[15,20]$.

From the perspective of optimization, and summarizing the geometrical
parameters of the prechamber in the tuple
$\mathbf{x}=(d_{\mathrm{bottle}},d_{\mathrm{bore}},h_{\mathrm{neck}})$, the
aim is to find a set of solutions to

$\max_{\mathbf{x}}k(\mathbf{x})\quad\text{s.t.}\quad|\vec{v}(\mathbf{x})|\leq
25.$ (1)

In BO [Fra18], both the objective $k$ and the constraint $|\vec{v}|$ are
considered as black box functions that are expensive to evaluate, and that are
thus approximated by surrogates. The most common surrogate models are Gaussian
processes [RW06]. Given a dataset
$\mathcal{D}=\\{(\mathbf{x}_{i},k(\mathbf{x}_{i}),|\vec{v}(\mathbf{x}_{i})|)\\}$
of different simulated settings, the Gaussian process models yield (Gaussian)
distributions $\hat{k}$ and $|\hat{\vec{v}}|$ for the target $k$ and the
constraint $|\vec{v}|$ for every possible parameterization $\mathbf{x}$. If
$\mathbf{x}$ is close to a point $\mathbf{x}_{i}$ in the dataset
$\mathcal{D}$, then the distribution of, e.g., $\hat{k}(\mathbf{x})$ has a
mean close to $k(\mathbf{x}_{i})$ and a small variance, whereas if
$\mathbf{x}$ is far from every point in $\mathcal{D}$, the variance of
$\hat{k}(\mathbf{x})$ will be large — there is high _epistemic uncertainty_.
During optimization, BO relies on these probabilistic models
$\hat{k}(\mathbf{x})$ and $|\hat{\vec{v}}(\mathbf{x})|$ to determine (using
acquisition functions, see [Fra18]) candidate geometries
$\mathbf{x}^{\bullet}$ that improve upon all geometries in $\mathcal{D}$ in
the sense of (1). The thus proposed candidate geometries are then evaluated
using CFD simulations, and consequently yield a new tuple
$(\mathbf{x}^{\bullet},k(\mathbf{x}^{\bullet}),|\vec{v}(\mathbf{x}^{\bullet})|)$.
This tuple is included in $\mathcal{D}$, and the probabilistic models
$\hat{k}(\mathbf{x})$ and $|\hat{\vec{v}}(\mathbf{x})|$ are updated. This
optimization loop is iterated until either a predetermined number of
iterations is exceeded or until the utility of the solution is satisfactory.

To address the optimization problem (1), we select $\mathbf{x}^{\bullet}$ via
the expected constrained improvement [GKX+14]. Specifically, if
$\mathbf{x}_{j}$ is the geometry in $\mathcal{D}$ that maximizes
$\hat{k}(\mathbf{x})=k(\mathbf{x})$, then $\mathbf{x}^{\bullet}$ is chosen as
an optimizer of

$\max_{\mathbf{x}}\mathbb{E}\left(\mathbf{1}(|\hat{\vec{v}}(\mathbf{x})|\leq
25)\cdot\max\\{0,\hat{k}(\mathbf{x})-k(\mathbf{x}_{j})\\}\right)\\\
=\max_{\mathbf{x}}PF(\mathbf{x})\cdot EI(\mathbf{x})$ (2)

where first factor is the probability that the candidate $\mathbf{x}$ is
feasible (i.e., that $|\hat{\vec{v}}(\mathbf{x})|\leq 25$) and where the
second factor is the standard expected improvement (w.r.t. $k$). Furthermore,
to exploit the fact that CFD simulations can be run in parallel, batch (or
parallel) BO was used, which responds with a set of $q$ candidate solutions
$\\{\mathbf{x}^{\bullet}_{1},\dots,\mathbf{x}^{\bullet}_{q}\\}$, cf. [Fra18,
eq. (14)]

In this work, BO was implemented using BoTorch [BKJ+20], a Python library that
contains functionality for constrained, multi-objective, and batch BO.
Specifically, both $\hat{k}$ and $|\hat{\vec{v}}|$ are modeled as Gaussian
processes with a fixed noise standard deviation of $0.005$ and a Matern kernel
with optimized hyperparameters. For optimal functionality of the BO, the three
geometry parameters $\mathbf{x}$ are scaled to the closed unit interval
$[0,1]$ and the objectives are scaled to unit variance ($k$ and $|\vec{v}|$)
and zero mean (only $k$). The expected constrained improvement was evaluated
using Monte Carlo sampling, with $1024$ samples. The batch size was set to
five.

A DOE study was performed with ten samples to create the initial data for
three consecutive BO iterations based on latin hypercube formulation. The
optimum value found by the DOE study was $k=160.38$ (see Table 1).

(a) DOE Initialization (10 data points)

(b) First BO Iteration (15 data points)

(c) Second BO Iteration (20 data points)

(d) Third BO Iteration (25 data points)

?figurename? 2: Progress of three consecutive BO iterations. The three columns
show the mean function of $\hat{k}(\mathbf{x})$ and its standard deviation for
two of the dimensions of $\mathbf{x}$ fixed to $\mathbf{x}^{\star}$,
respectively

##  3 Results

The results for the optimum values found by the DOE initialization and three
consecutive BO iterations can be found in Table 1. While the first BO
iteration shows considerable improvement ($k=246.46$) compared to the DOE
initialization, the second BO iteration was not able to further increase the
objective within the following five candidate positions ($k=244.56$). The
third and final BO iteration again yields a new optimum value ($k=263.16$)
while still fulfilling the constraint ($|\vec{v}|=22.53$).

Figure 2 illustrates the progress of BO over the three iterations. Given that
after the third iteration the optimal value in $\mathcal{D}$ is
$\mathbf{x}^{\star}=(d_{\mathrm{bottle}}^{\star},d_{\mathrm{bore}}^{\star},h_{\mathrm{neck}}^{\star})$,
the three columns show the mean function of $\hat{k}(\mathbf{x})$ and its
standard deviation for two of the dimensions of $\mathbf{x}$ fixed to
$\mathbf{x}^{\star}$, respectively. One can see that by including more
training data, the surrogate $\hat{k}$ becomes more and more structured and
that, especially w.r.t. $d_{\mathrm{bottle}}$, the epistemic uncertainty
reduces. The working of BO is best illustrated in the third and fourth rows of
the figure: Figure 2(c) depicts the surrogate model $\hat{k}$ trained on 20
data points. The vertical dashed lines indicate the candidate position
$\mathbf{x}^{\bullet}$ that maximizes the acquisition function (2). It can be
seen that the candidate position lies at a local optimum of the mean of
$\hat{k}$ for $d_{\mathrm{bore}}$, with a value close to 180. Running CFD
simulations with the thus suggested geometry parameters $\mathbf{x}^{\bullet}$
yields a value of $k(\mathbf{x}^{\bullet})$ exceeding 250. Figure 2(d) then
shows the surrogate $\hat{k}$ after including $\mathbf{x}^{\bullet}$ (and the
remaining four geometries of the batch) in the training data $\mathcal{D}$. It
can be seen that training was successful in the sense that the
$\mathbf{x}^{\bullet}=\mathbf{x}^{\star}$ now lies on the mean curve of
$\hat{k}$ and that the epistemic uncertainty at this position is small.

?tablename? 1: Optimum values found by the DOE initialization (10 data points)
and three consecutive BO iterations (5 data points per iteration).

| $k$ | $|\vec{v}|$ | $d_{\mathrm{bottle}}$ | $d_{\mathrm{bore}}$ | $h_{\mathrm{neck}}$
---|---|---|---|---|---
DoE Init. | 160.38 | 17.93 | 10.20 | 0.89 | 18.75
$1^{\mathrm{st}}$ BO It. | 246.46 | 22.70 | 10.02 | 0.88 | 18.26
$2^{\mathrm{nd}}$ BO It. | 244.56 | 22.11 | 9.90 | 0.88 | 18.18
$3^{\mathrm{rd}}$ BO It. |  263.16 |  22.53 |  11.81 |  0.88 |  19.74
|  |  |  |  |

Despite the fact that BO is a global optimization method, we have observed in
our experiments that most candidate geometries were suggested in the vicinity
of optima of the surrogate function. In other words, BO _exploited_ knowledge
about the currently best geometry rather than _exploring_ geometries that have
not been simulated so far. Upon inspection of Figure 2, one can see that the
target $k(\mathbf{x})$ varies substantially with the bore diameter; with only
few training data, this variation cannot be captured in the surrogate
$\hat{k}(\mathbf{x})$. (Indeed, the surrogate $\hat{k}$ in Figure 2(c)
underestimated the value of the $k$ at the candidate position indicated with
the dashed lines, as can be seen by the star in Figure 2(d).) This is
aggravated by the fact that the kernel parameters are optimized based on
$\mathcal{D}$ and, thus, may not be optimal for the black box functions $k$
and $|\vec{v}|$. Indeed, as the center image in Figure 2(c) suggests, the
automatic hyperparameter optimization may have terminated with a kernel width
that is too large to capture the variation of $k$ with the bore diameter
(especially when compared with the corresponding images in Figure 2(b) and
2(d), respectively). A manual setting of these hyperparameters, however, is
not possible since the functions $k$ and $|\vec{v}|$ are unknown a priori.

Figure 3 shows the determined $k$-field in the cross section of the prechamber
which allows a detailed analysis of the BO design parameter optimization
results. While Figure 3a and Figure 3b represent an unfavourable and a
favourable solution of the initial DOE set, respectively, Figure 3c shows the
final resulting parameter combination found by the BO. It can be clearly
observed, that the formation of a vortex in the prechamber is promoted by the
appropriate combination of the chosen design parameters. Since the mass
flowing from the main chamber into the prechamber is mainly influenced by the
piston movement, a smaller value of $d_{\mathrm{bore}}$ results in higher
overflow velocity and consequently higher turbulence velocities in the
prechamber for the predefined parameter range. It should be mentioned at this
point that this effect is not unlimited since the reduction of the overflow
area can lead to choking effects. The combination a small $d_{\mathrm{bore}}$
value and high values of $d_{\mathrm{bottle}}$ as well as $h_{\mathrm{neck}}$
leads to a well-developed vortex with maximum values of $k$ in the vicinity of
the spark plug.

(a)(b)(c)

?figurename? 3: Turbulence intensity $k$-field in the prechamber cross section
of one unfavourable design combination from the initial set (a), the best
design combination of the initial set (b) and the final BO predicted design

##  4 Summary and Outlook

In the present study, a combination of CFD and BO was applied to optimize the
design of a large bore gas engine prechamber. Three design parameter namely
the diameter of the prechamber bottle, the diameter of the overflow bores and
the height of the prechamber neck were chosen as input parameters. The target
of the optimization process was maximizing the turbulence intensity $k$ under
a velocity magnitude constraint in the vicinity of the spark plug. After the
initialization of the system via ten DOE-based samples, three BO steps were
carried out. The analysis of the BO indicate that both the strong variability
of $k$ with the bore diameter and the effects of automatic hyperparameter
optimization appear to slow down BO. One may therefore assume that
substantially more iterations are required for BO to enter the exploration
phase and, consequently, suggest geometry parameters distinct from the current
optima. An option to speed up BO is to chose acquisition functions that permit
influencing the trade-off between exploration and exploitation, such as the
upper confidence bound. Future work shall be devoted to investigate these
alternatives. Furthermore, additional investigations regarding the complexity
of the physical system and its sensitivity to design parameter changes will be
performed in order to increase the input and output parameter space.

## Acknowledgments

The authors acknowledge the financial support of the Austrian COMET —
Competence Centers for Excellent Technologies — Programme of the Austrian
Federal Ministry for Climate Action, Environment, Energy, Mobility, Innovation
and Technology, the Austrian Federal Ministry for Digital and Economic
Affairs, and the States of Styria, Upper Austria, Tyrol, and Vienna for the
COMET Centers LEC EvoLET and Know-Center, respectively. The COMET Programme is
managed by the Austrian Research Promotion Agency (FFG).

##  ?refname?

  * [BKJ+20] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy.  BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization.  In Advances in Neural Information Processing Systems 33, 2020.
  * [BNGS+20] J. Benajes, R. Novella, J. Gomez-Soriano, I. Barbery, C. Libert, F. Rampanarivo, and M. Dabiri.  Computational assessment towards understanding the energy conversion and combustion process of lean mixtures in passive pre-chamber ignited engines.  Applied Thermal Engineering, 178:115501, 2020.
  * [Fra18] Peter I. Frazier.  A Tutorial on Bayesian Optimization, 2018.  arXiv:1807.02811.
  * [FZQ+18] Liyan Feng, Jun Zhai, Chuang Qu, Bo Li, Jiangping Tian, Lei Chen, Weiyao Wang, Wuqiang Long, and Bin Tang.  The influence of the enrichment injection angle on the performance of a pre-chamber spark ignition natural-gas engine.  Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering, 232(5):679–694, 2018.
  * [GBV+19] Elisa Guelpa, Aldo Bischi, Vittorio Verda, Michael Chertkov, and Henrik Lund.  Towards future infrastructures for sustainable multi-energy systems: A review.  Energy, 184:2–21, 2019.  Shaping research in gas-, heat- and electric- energy infrastructures.
  * [GBY+21] Haiwen Ge, Ahmad Hadi Bakir, Suren Yadav, Yunseon Kang, Siva Parameswaran, and Peng Zhao.  Cfd optimization of the pre-chamber geometry for a gasoline spark ignition engine.  Frontiers in Mechanical Engineering, 6:599752, 2021.
  * [GG21] Clemens Gößnitzer and Shawn Givler.  A new method to determine the impact of individual field quantities on cycle-to-cycle variations in a spark-ignited gas engine.  Energies, 14(14):4136, 2021.
  * [GKK+19] M. Ghaderi Masouleh, K. Keskinen, O. Kaario, H. Kahila, S. Karimkashi, and V. Vuorinen.  Modeling cycle-to-cycle variations in spark ignited combustion engines by scale-resolving simulations for different engine speeds.  Applied Energy, 250:801–820, 2019.
  * [GKX+14] Jacob R Gardner, Matt J Kusner, Zhixiang Eddie Xu, Kilian Q Weinberger, and John P Cunningham.  Bayesian optimization with inequality constraints.  In Proc. Int. Conf. on Machine Learning (ICML), pages 937–945, 2014\.
  * [GPP22] Clemens Gößnitzer, Stefan Posch, and Gerhard Pirker.  Development of a framework for internal combustion engine simulations in openfoam.  In 8th European Congress on Computational Methods in Applied Sciences and Engineering: ECCOMAS CONGRESS 2022, 2022.
  * [NPGS+20] Ricardo Novella, Jose Pastor, Josep Gomez-Soriano, Ibrahim Barbery, Cedric Libert, Fano Rampanarivo, Chistou Panagiotis, and Maziar Dabiri.  Experimental and numerical analysis of passive pre-chamber ignition with egr and air dilution for future generation passenger car engines.  In WCX SAE World Congress Experience. SAE International, apr 2020\.
  * [Opea] OpenFOAM Development Repository.  Github commit c95d964 from 20th November 2022, https://github.com/OpenFOAM/OpenFOAM-dev.  Accessed: 2023-01-03.
  * [Opeb] The OpenFOAM Foundation.  https://openfoam.org/.  Accessed: 2023-01-03.
  * [PGO+22] Stefan Posch, Clemens Gößnitzer, Andreas B Ofner, Gerhard Pirker, and Andreas Wimmer.  Modeling cycle-to-cycle variations of a spark-ignited gas engine using artificial flow fields generated by a variational autoencoder.  Energies, 15(7):2325, 2022.
  * [PWZ+21] S. Posch, H. Winter, J. Zelenka, G. Pirker, and A. Wimmer.  Development of a tool for the preliminary design of large engine prechambers using machine learning approaches.  Applied Thermal Engineering, 191:116774, 2021.
  * [RF03] R.P. Roethlisberger and D. Favrat.  Investigation of the prechamber geometrical configuration of a natural gas spark ignition engine for cogeneration: part ii. experimentation.  International Journal of Thermal Sciences, 42(3):239–253, 2003\.
  * [RW06] Carl Edward Rasmussen and Christopher K.I. Williams.  Gaussian Processes for Machine Learning.  The MIT Press, 2006.
  * [SMB+22] Mickael Silva, Balaji Mohan, Jihad Badra, Anqi Zhang, Ponnya Hlaing, Emre Cenker, Abdullah S. AlRamadan, and Hong G Im.  Doe-ml guided optimization of an active pre-chamber geometry using cfd.  International Journal of Engine Research, 0(0):14680874221135278, 2022.
  * [SSH+20] Mickael Silva, Sangeeth Sanal, Ponnya Hlaing, Bengt Johansson, Hong G. Im, and Emre Cenker.  Effects of geometry on passive pre-chamber combustion characteristics.  In WCX SAE World Congress Experience. SAE International, apr 2020\.
  * [WSP+20] Hubert Winter, Eduard Schneßl, Gerhard Pirker, Jan Zelenka, and Andreas Wimmer.  Application of CFD simulation to optimize combustion in prechamber gas engines with port injection.  In Proceedings of the 6th European Conference on Computational Mechanics: Solids, Structures and Coupled Problems, ECCM 2018 and 7th European Conference on Computational Fluid Dynamics, ECFD 2018, pages 894–905, 2020.
  * [ZALY22] Sipeng Zhu, Sam Akehurst, Andrew Lewis, and Hao Yuan.  A review of the pre-chamber ignition system applied on future low-carbon spark ignition engines.  Renewable and Sustainable Energy Reviews, 154:111872, 2022.
  * [ZYE+20] Anqi Zhang, Xin Yu, Nayan Engineer, Yu Zhang, and Yuanjiang Pei.  Numerical investigation of pre-chamber jet combustion in a light-duty gasoline engine.  In Internal Combustion Engine Division Fall Technical Conference, volume 84034, page V001T06A013. American Society of Mechanical Engineers, 2020.