#  Neural network analysis of neutron and X-ray reflectivity data: automated
analysis using _mlreflect_ , experimental errors and feature engineering

Alessandro Grecoa Vladimir Starostina Evelyn Edela Valentin Munteanua Nadine
Rußeggera Ingrid Daxa Chen Shenb Florian Bertramb Alexander Hinderhofera
Alexander Gerlacha Frank Schreibera [a] Institute of Applied Physics, Auf der
Morgenstelle 10, University of Tübingen, 72076 Tübingen, Germany
Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

###### Abstract

This work demonstrates the Python package _mlreflect_ which implements an
optimized pipeline for the automized analysis of reflectometry data using
machine learning. The package combines several training and data treatment
techniques discussed in previous publications. The predictions made by the
neural network are accurate and robust enough to serve as good starting
parameters for an optional subsequent least mean squares (LMS) fit of the
data. It is shown that for a large dataset of 242 reflectivity curves of
various thin films on silicon substrates, the pipeline reliably finds an LMS
minimum very close to a fit produced by a human researcher with the
application of physical knowledge and carefully chosen boundary conditions.

Furthermore, the differences between simulated and experimental data and their
implications for the training and performance of neural networks are
discussed. The experimental test set is used to determine the optimal noise
level during training. Furthermore, the extremely fast prediction times of the
neural network are leveraged to compensate for systematic errors by sampling
slight variations of the data.

###  I Introduction

X-ray and neutron reflectometry (XRR and NR) are established surface
scattering techniques that are routinely used to characterize solid and liquid
thin films [1, 2, 3, 4]. They offer a non-invasive way to determine the
structural, morphological or magnetic properties of a large variety of samples
[5, 6, 7] and can also be employed in real-time for in situ measurements [8].
For decades, the conventional way to analyze reflectivity data has been the
iterative least mean squares (LMS) or $\chi^{2}$ fitting of the data with a
theoretical model [9, 10, 11]. However, due to the well-known phase problem in
scattering, the reconstruction of the scattering length density (SLD) profile
from the reflectivity data is inherently ambiguous. This means that this
method typically requires significant expertise and prior knowledge about the
system, since for all but the simplest cases, there exist many possible
solutions. Furthermore, even when the solution space is restricted, finding
the global minimum is usually very time-consuming due to several local minima
on the mean squared error (MSE) surface. For this purpose, various software
packages have been developed over the years that use sophisticated
minimization algorithms [12, 13, 14, 15, 16, 17]. However, all of these
approaches are iterative in nature and thus, usually computationally slow.
Recently, machine learning-based methods have been proposed, that could avoid
a lengthy search of the MSE surface, by providing an immediate guess for the
thin film parameters that is already very close to the ground truth [18, 19,
20, 21, 22] or by encoding the reflectometry data into a latent space where
the error surface does not have as many local minima [23].

This paper demonstrates a Python-based reflectivity data analysis pipeline
called _mlreflect_ that combines a fully-connected neural network regressor
with several preprocessing and postprocessing steps to reliably predict the
thickness, roughness and SLD of a thin film layer. While the principle of the
neural network itself and the preprocessing have been discussed previously
[18, 22], here we focus on the differences between simulated and experimental
data and show how this knowledge can be used to further optimize the obtained
results. We tested the performance of the pipeline on a large experimental
dataset of 242 XRR curves from different samples by comparing the result of
the pipeline with manually supervised LMS fits that include physical knowledge
and carefully chosen boundary conditions. This is a quantitative and
qualitative difference compared to other similar studies, where most of the
performance analysis is done with simulated data. In this context, we discuss
the effect that experimental deviations from the theory can have on the
training and prediction quality of the neural network. Using an example curve,
we show how the extremely fast prediction speed of the neural network can also
be leveraged to compensate for small experimental errors.

###  II Description of the analysis pipeline

Figure 1: A schematic description of the analysis pipeline. The pipeline
consists of three main steps: I. preprocessing, II. parameter prediction via
the neural network and III. postprocessing. Step I. includes geometrical and
other experiment-specific corrections. The data is also normalized,
transformed into $q_{z}$ space, interpolated and standardized. In step II.,
the preprocessed data is fed into a trained, fully-connected neural network
that yields an initial guess for the thin film parameters. During step III.,
this initial guess is used as starting parameters for a fast Levenberg-
Marquardt fit that finds the nearest LMS minimum.

Our proposed analysis pipeline _mlreflect_ is fully written in Python and is
available as open source on GitHub (https://github.com/schreiber-
lab/mlreflect) and can also be downloaded directly from the Python Package
Index (https://pypi.org/project/mlreflect/). The Supporting Information to
this manuscript contains a step-by-step tutorial in the form of executable
Jupyter notebooks (also available as PDF version). In addition, the tutorial,
installation instructions and a full API documentation of the _mlreflect_
package are hosted on https://mlreflect.readthedocs.io/en/latest/. The neural
network itself is implemented using TensorFlow [24]. It uses the matrix
formalism implemented in the refl1d package [13] to simulate the reflectivity
data. The workflow of the package can be conceptually separated into three
steps: I. preprocessing, II. prediction and III. postprocessing, as depicted
in Figure 1. Each of these steps is described in the following.

During step I., the reflectivity data is automatically read from its raw
format and several types of preprocessing procedures are applied. First, the
raw data is converted into the standard $R(q_{z})$ format where $R$ is the
normalized reflected intensity and $q_{z}$ the momentum transfer vector along
the surface normal. The preprocessing operations necessary depend on how the
raw data is saved, but usually the data has to be corrected in some form. In
our case, the raw data contains the reflected intensity at different
scattering angles that must be corrected for the varying beam attenuation at
different angles. Then, the intensity is corrected to account for the changing
beam footprint on the sample at different angles, which amounts to a
multiplication of the data with a geometric factor [25]. Here we assume a flat
sample and a beam with a Gaussian profile but, in principle, corrections for
other sample or beam shapes can be implemented at this stage. The data is then
normalized by dividing by the highest intensity value and transformed from
angular space into $q_{z}$ space.

After that, the intensity values are interpolated on a logarithmic scale to
109 equally-spaced $q_{z}$ points ranging from 0.02 to
$0.15\text{\,}{\mathrm{\text{Å}}}^{-1}$, which corresponds to the input size
of the neural network. Lastly, the data is standardized in the same way as was
described in [22], which ensures that each value of the input vector is on a
similar scale. The effect on the general shape of the curves is comparable to
multiplying the data with the inverse of the Fresnel reflectivity
$R_{\text{F}}(q_{z})\propto q_{z}^{-4}$, but, importantly, avoids the
divergence for small values of $q_{z}$, i.e. close to and below the total
reflection edge (TRE), where the kinematical approximation does not hold [26].

To obtain the initial parameter prediction (step II. in Figure 1) the
preprocessed input vector is fed into the trained neural network model. The
neural network is a fully-connected model that takes an input of 109 discrete
intensity points and outputs 3 thin film parameters: the film thickness, the
Névot-Croce film roughness [27] and the real part of the scattering length
density of the film. The model has 3 hidden layers with 512 neurons each. The
training loss was calculated as the mean squared error between the normalized
predicted and ground truth parameters. This architecture is similar to what
has been described in the literature before [18, 22, 20], but to reduce
training and inference times the number of parameters was reduced. The model
was trained with 250000 simulated reflectivity curves with a batch size of
512. For every batch, uniform noise and curve scaling were applied to each
curve to avoid overfitting as described before [22]. The optimal noise level
during training was identified to be 0.3, which will be discussed in more
detail later. Finally, the inputs were standardized as described above.

The training data was generated assuming a sample structure consisting of a
thin film on top of an oxide-capped silicon substrate with air as an ambient
medium and with X-rays as the probe. The thin film parameters in the training
data spanned a large range of 20–$1000\text{\,}\mathrm{\text{Å}}$ for the
thickness, 0–$100\text{\,}\mathrm{\text{Å}}$ for the roughness and
1–$14\text{\times}{10}^{-6}\text{\,}{\mathrm{\text{Å}}}^{-2}$ for the SLD.
Furthermore, we restricted the roughness to values no higher than half the
thickness since these scenarios are not well described by the theoretical
model used. We note that a similar approach could easily be employed for
neutrons or other sample structures by retraining the neural network with
different training data. We also expect this approach to work for a larger
number of layers as long as the trained parameter space does not create too
many ambiguous solutions, i.e., the number and range of fitting parameters
should remain similar. For a larger parameter space, a larger $q_{z}$ range
might be necessary to reduce ambiguity in the data. In our case, the $q_{z}$
range was limited to avoid conflicts with the Bragg peaks of organic molecules
around $0.3\text{\,}{\mathrm{\text{Å}}}^{-1}$ which are not described by the
slab model.

Lastly, during step III., the initially predicted thin film parameters are fed
into an LMS minimizer to obtain the parameters which produce the best fit.
Since the initial predictions are already very close to the ground truth, we
chose a simple Levenberg-Marquardt minimizer [28] over a more powerful, but
slower algorithm.

###  III Performance test on thin films

The performance of the analysis pipeline was tested on 242 experimental XRR
curves from in situ and ex situ experiments of 9 organic thin films on Si/SiOx
(1–79 curves per sample at different thicknesses). The distributions of
thickness, roughness and SLD of the film within this test set is shown in
Figure 2. The measurements were conducted using three different synchrotron
radiation sources, i.e. the ESRF [29], DESY [30] and the SLS [31], as well as
using our own laboratory source. To obtain a benchmark, each reflectivity
curve was first fitted on a logarithmic scale with an LMS fit based on the
commonly used differential evolution algorithm [32] using manually chosen
initial values and bounds for each parameter. The thin film model used for the
fit was the same as what was used for training the neural network. In the
following analysis, we assume that these manually fitted parameters represent
the “ground truth” and thus the performance of our pipeline will be measured
as the absolute error with respect to this ground truth.

Figure 2: Ground truth distribution of the three sample parameters thickness
(a), roughness (b) and SLD (c) within the experimental test set of 242 XRR
curves. The parameters were obtained by a conventional LMS fit.

In the following, we compare the ground truth with prediction results of the
neural network as well as the results of an automized subsequent LMS fit using
the predicted parameters. Across all 242 curves, the neural network
predictions have a median absolute error (median percentage error) of
$6.0\text{\,}\mathrm{\text{Å}}$ (7.1%) for the film thickness,
$2.0\text{\,}\mathrm{\text{Å}}$ (12.4%) for the interface roughness and
$0.72\text{\times}{10}^{-6}\text{\,}{\mathrm{\text{Å}}}^{-2}$ (6.8%) for the
SLD. This is a significant improvement to our first published model [18], both
on an absolute scale as well as on a relative scale, since the possible ranges
for the thickness and roughness parameter have been greatly expanded. Thus,
the network is generalized over a larger parameter space compared to
previously published results. We note that since all of our data stems from
organic thin films, the SLDs in the test set are mainly clustered around
10–$13\text{\times}{10}^{-6}\text{\,}{\mathrm{\text{Å}}}^{-2}$. Nevertheless,
we assume that our results are not specific to the SLD range of the test data,
since the network was trained equally with SLDs ranging from
1–$14\text{\times}{10}^{-6}\text{\,}{\mathrm{\text{Å}}}^{-2}$. We also want to
highlight that the dataset also contains curves with a high roughness-to-
thickness ratio where the Kiessig oscillations are strongly damped. Among the
emerging solutions offered in this field, discussions about the performance on
curves with little to no features are mostly absent. This is of course due to
the challenge of extracting information from data that inherently contains
less information. Yet, the network presented here also performs well on
experimental data with high relative roughness.

Figure 3: a) Box plot of the absolute errors for 242 measured reflectivity
curves for each of the three predicted parameters. The upper and lower edge of
the boxes represent the 1st and 3rd quartile with the horizontal line inside
the boxes denoting the median. The blue boxes represent the error compared to
the pure neural network predictions. The red boxes represent the error after
applying a simple LMS minimization using the neural network predictions as
starting parameters. The green boxes show the error for the case when a
$q_{z}$ shift optimization has been performed before the LMS fit. b) The same
box plots of the median error but as a percentage of the ground truth.
All results were obtained for a training noise level of $n=0.3$.

The next step in the pipeline is to further refine these results via an LMS
fit using the predictions from the neural network as starting parameters.
Since the predictions are robust and already quite close to the ground truth
there is no need for powerful but slow minimization algorithms such as genetic
or differential evolution algorithms, which are normally employed to find the
global minimum. Thus, finding the minimum takes only a fraction of a second
per curve and can be fully automized. After this refinement procedure, the
median absolute error (median percentage error) was even closer to the ground
truth with $2.3\text{\,}\mathrm{\text{Å}}$ (2.3%) for the thickness,
$1.0\text{\,}\mathrm{\text{Å}}$ (5.8%) for the roughness and
$0.47\text{\times}{10}^{-6}\text{\,}{\mathrm{\text{Å}}}^{-2}$ (4.3%) for the
SLD. A comparison of the error distributions before and after refinement is
shown in Figure 3. A detailed breakdown of the prediction error with respect
to each parameter can be found in Figures S2–S10 of the Supporting
Information.

The residual error can be attributed to the fact that every fit has a finite
accuracy and hence the ground truth itself contains a certain error. We
roughly estimate this error to be at least $\pm 10\,\%$ for each parameter,
which would be comparable to the reported error of the neural network. Thus,
these results show that the analysis pipeline as described above performs
similarly to a human researcher in most circumstances. Furthermore, it is
important to note that the results were obtained much faster than via a human-
guided fit. Excluding the time it took to train the neural network (about 20
minutes for a given sample structure), the initial parameter predictions of
the 242 curves were obtained after only 1 second with about 2 additional
minutes for the further refinement steps, resulting in a total fitting time of
about 0.4s per curve. In contrast, producing the ground truth fits took about
6 hours because of the need to carefully select fitting boundaries to prevent
the fit from running into non-physical minima.

###  IV Differences between simulated and experimental data

A well-known property of artificial neural networks is that they require large
amounts of representative training data to learn a generalized model and not
overfit the training set. In the context of the work presented here, i.e.
supervised learning using scattering data, this would mean acquiring thousands
of scattering patterns from a large variety of different samples and analyzing
them manually to create the training set. Since this is a quite time-consuming
and challenging task, neural network models in the field of scattering physics
are typically trained with simulated data based on well-established
theoretical models. In most cases, the simulation is additionally modified
with certain artifacts, such as noise, to better mimic experimental
conditions. However, to what degree this is necessary is difficult to estimate
since the only available metric is typically the performance on other
simulated data (validation loss), which is expected to decrease with
increasing perturbations.

Figure 4: Comparison of the loss calculated from a simulated test set (100000
curves) and an experimental test set (242 curves) for different levels of
uniform noise that were applied to the training data. For each noise level a
separate model was trained. With increasing noise level, the loss from the
simulated data increases linearly while the loss from the experimental data
shows a clear minimum at noise levels 0.3–0.35.

In this study, we investigated how applying uniform noise to the training data
affects the neural network performance on our large experimental dataset of
242 curves. We trained 11 copies of the same neural network (as described
above) with training data with different noise levels $n$ where each data
point $R_{i}^{*}$ in the noisy curve was sampled uniformly between the values
$R_{i}(1-n)$ and $R_{i}(1+n)$. Thus, $n$ denotes the maximum relative change
of a given data point $R_{i}$ of a given simulated curve. The $n$ for each
trained model ranged from 0 to 0.5 in 0.05 increments. It is important to note
that the applied uniform noise is not meant to model a specific physical noise
type, such as Poisson noise for counting statistics. Rather, uniform noise was
chosen as a $q$-independent catch-all noise that affects the whole curve
equally and thus, makes the neural network robust against errors across the
entire $q$ range.

Figure 4 shows a comparison of the losses calculated with a simulated test set
as well as with the experimental test set for each model. Since the loss is
calculated as the mean squared error of all three (normalized) sample
parameters, it is a unitless measure for the accuracy of the model. For $n=0$,
the simulated test set shows a loss close to zero ($\sim 10^{-7}$), whereas
the loss based on the experimental data is about 5 orders of magnitude higher.
This shows that without any noise, the neural network significantly overfits
the simulation and thus performs suboptimally on real data. As expected, the
loss of the simulated data increases monotonically with increasing noise.
However, the performance on the real data improves significantly with
increasing noise up to a noise level of 0.3–0.35. Beyond this, even higher
noise levels seem to again worsen the performance. This very clearly
demonstrates that there exists an optimal noise level for which the added
noise acts as an effective regularization technique that prevents overfitting.
If the noise level is too high, however, the consequent lack of information is
likely detrimental to the training. Thus, we identified $n=0.3$ to be the
ideal noise level for data similar to our testing set, which notably contains
data from different X-ray sources. Furthermore, Figure S1 of the Supporting
Information shows that the optimal training noise does not significantly
change for subsets with different noise levels (0.1–0.5) within the
experimental test set. Thus, we set the default value of the noise level in
our analysis pipeline to 0.3. Datasets that differs significantly from our
test set in terms of experimental artifacts might of course produce slightly
different results, although we expect the general trend to be the same. This
highlights the importance of having a large experimental test set with
representative experimental artifacts since metrics only based on simulated
data are clearly not sufficient to evaluate the training progress.

###  V Influence of systematic measurement errors

All reflectometry measurements are performed with a finite accuracy due to
various error sources. These errors are detrimental to the experiment and can
impede the extraction of information from the data and therefore should be
avoided or minimized as much as possible. However, a finite error inevitably
remains for every measurement. Among the possible statistical errors are the
signal-to-noise ratio, the angular resolution of the diffractometer and the
spectral resolution of the source. Among the systematic errors are, for
example, the convolution of the data with the slit functions, the accuracy of
the sample alignment and the accuracy of the footprint correction (i.e. how
accurate the beam and sample shape can practically be determined).

Having imperfect data obviously impacts the analysis, since the data deviates
from the ideal physical model it is compared with. Since the neural network
model presented here is trained to solve a very particular task that assumes
well-defined data, these errors can negatively impact the prediction quality.
In general, it is easier to make the neural network robust against statistical
errors by introducing them during training, as described before. However,
sometimes systematic errors, such as a small misalignment can also seriously
misguide the ML prediction, as shown in Figure 5. Therefore, it would be
useful to correct or compensate some of these errors during inference time
after the data has been acquired.

As a solution, we propose an automated method for sampling through slight
variations of the data, exploiting both the sensitivity and speed of our
neural network model. Since the neural network assumes data that conforms to
an idealized physical model, it might fail if the data contains anomalies with
respect to that model. Since predictions with the neural network are very
fast, it is possible to scan through thousands of modified reflectivity curves
within less than a second. For each of these variants, the log MSE between the
data and the predicted curve can be calculated and only the one with the
lowest error is subsequently selected. We demonstrate an implementation of
this method that identifies small systematic alignment errors and
automatically applies an appropriate shift to the data.

Figure 5a shows an XRR measurement of a $690\text{\,}\mathrm{\text{Å}}$ thick
N,N’-Dioctyl-3,4,9,10-perylenedicarboximide (PDI-C8) film on Si/SiOx which was
measured and tested in addition to the 242 test curves. Here, in contrast to
the previously shown test set, the normal pipeline as described above did not
converge to the correct minimum. The reason for this is the much higher
thickness of the film, which leads to denser Kiessig oscillations in the data.
This, in turn, creates many narrow minima on the MSE surface for the LMS
algorithm to get trapped in. As a result, the neural network prediction needs
to be even closer to the ground truth for the subsequent fit to converge.
Table 1 shows the predicted thin film parameters in comparison to the ground
truth. A possible reason for the suboptimal neural network prediction might be
small imperfections in the data due to finite measurement errors, such as a
small variation in sample alignment. In regions of high derivatives, even a
small shift of the data along the $q_{z}$ axis can lead to strong differences
in the observed intensities at a given $q_{z}$ value, even on a logarithmic
scale. Of course, if the data has dense oscillations, this effect becomes more
pronounced. For models trained on simulated data, this can be critical, since
normally a substantial change of certain input neurons, especially near the
TRE, corresponds to important information and will be interpreted by the
network accordingly. To check whether this can be remedied, we shifted the
$q_{z}$ values during the interpolation step by a small value $\Delta q_{z}$
and repeated the prediction. This was done 1000 times with randomly sampled
$\Delta q_{z}$ ranging from $-1\text{\times}{10}^{-3}\text{\,}$ to
$1\text{\times}{10}^{-3}\text{\,}{\mathrm{\text{Å}}}^{-1}$. Then, for each
prediction, the quality of the prediction was evaluated by calculating the log
MSE between the corresponding simulation and the measured curve.

Figure 5: a) Comparison of the neural network predictions from reflectivity data from a $690\text{\,}\mathrm{\text{Å}}$ thick PDI-C8 film on Si/SiOx. The blue curve shows the native prediction whereas the red curve shows the prediction after the data was shifted by $\Delta q_{\text{min}}=$5.2\text{\times}{10}^{-4}\text{\,}{\mathrm{\text{Å}}}^{-1}$$ before the interpolation step. It is apparent that the latter is in much better agreement with the data. b) Shows the log MSE between the predicted curve and the data for different $\Delta q_{z}$. The minimum MSE at $\Delta q_{\text{min}}$ is indicated by the dashed line. Table 1: Predicted and fitted thin film parameters based on the reflectivity data of a PDI-C8 film on Si/SiOx (shown in Figure 5). The ground truth labels were obtained via a manually supervised LMS fit. After applying the described $q_{z}$ variation, the prediction results improved significantly. A subsequent LMS refinement only led to comparatively small improvements. | Thickness [Å] | Roughness [Å] | SLD [$10^{-6}$${\mathrm{\text{Å}}}^{-2}$]
---|---|---|---
ground truth | 688.3 | 27.1 | 10.5
prediction | 536.7 | 30.3 | 11.2
shift + prediction | 690.8 | 31.0 | 11.0
shift + prediction + fit | 690.5 | 27.5 | 10.8

When plotting the log MSE between the prediction and the input against $\Delta
q_{z}$ (Figure 5), we observed a value $\Delta
q_{\text{min}}=$5.2\text{\times}{10}^{-4}\text{\,}{\mathrm{\text{Å}}}^{-1}$$
for which the log MSE shows a clear minimum. From Figure 5a it is apparent
that the predicted curve based on the shifted data shows much better agreement
with the data than the normal prediction. The corresponding predicted
parameters for $\Delta q_{\text{min}}$ are shown in Table 1 and are much
closer to the ground truth values (comparable to values shown in the previous
section). This indicates that there exists a certain shift $\Delta
q_{\text{min}}$ that can (at least partially) compensate for the experimental
error. This is especially valuable since it allows the pipeline to continue
with the LMS refinement step, which ultimately leads to a near-perfect fit.

It is interesting to note that $\Delta q_{\text{min}}$ is very small,
corresponding to a change of the angle of incidence of only about
$4\text{\times}{10}^{-3}\text{\,}$ degrees for a wavelength of
$1.54\text{\,}\mathrm{\text{Å}}$. It seems intuitive that such a small shift
in the data could be caused by a variety of the above mentioned error sources.
However, although $\Delta q_{\text{min}}$ is seemingly small, due to the high
derivatives close to the TRE and the Kiessig fringes, shifting the data by
$\Delta q_{\text{min}}$ still has a noticeable effect on each data point. For
conventional LMS fitting, this might not seem critical at first, since the MSE
surface likely has a minimum close to the real one in terms of the film
thickness. However, for the roughness and density parameters this might not be
the case and thus, most fitting programs allow the user to manually shift the
data if necessary.

While in principle any type of modification like this could be conceivably
applied to the data to scan for the lowest MSE, we observed significantly
better results with this method rather than, for example, adding Gaussian
noise. This is because a translation of the curve preserves most of the
information in the data while still varying every data point, in contrast to
Gaussian sampling which is $q$-independent and inevitably destroys
information.

To test the stability of this method, we applied the $\Delta q_{z}$ sampling
procedure to all 242 curves discussed in the previous section (where the
pipeline already succeed) and compared the results with the original mean
absolute error. When looking at Figure 3, it becomes clear that scanning for
$\Delta q_{\text{min}}$ did not harm the mean absolute error, but instead even
improved the results slightly for all three parameters. While the log MSE of
the predictions is already very close to the minimum, most of the data likely
still has a finite alignment error which, however, was not sufficient to
affect the prediction. Hence, this could still be compensated by applying a
small shift, ultimately leading to an even better fit. Because this screening
for $\Delta q_{z}$ yielded significant improvements on some data and was
relatively fast, we decided to routinely add this to the analysis pipeline.

###  VI Fourier transforms as a method for feature engineering

The specular reflectivity from a single layer on a substrate well above the
critical angle can be approximately described by

$R(q_{z})=R_{\text{F}}(q_{z})\left|\int_{-\infty}^{\infty}\frac{d\rho(z)}{dz}e^{iq_{z}z}dz\right|^{2}$
(1)

i.e. the product of the Fresnel reflectivity from a flat surface and the
squared Fourier transform of the SLD contrast of the sample along the surface
normal [26, 33]. Although the phase of the Fourier transform is lost by taking
its absolute square, the inverse Fourier transform of
$R(q_{z})/R_{\text{F}}(q_{z})$ still carries some important information such
as the frequency of the Kiessig oscillations (and thus the film thickness). As
a result, performing an inverse Fourier transform on the reflectivity data
presents itself as an obvious way to create additional input features that may
facilitate the neural network training.

To test this hypothesis, we trained a neural network model with an additional
preprocessing step before the first layer that performs a Fast Fourier
transform on the standardized input and adds the real and imaginary Fourier
components to it, leading to a input layer size of 219 neurons. All other
model parameters and training ranges were kept the same as described above.
When testing the trained model on the 242 experimental curves, we found that
the model performed similarly to the model without the added Fourier
transform. The median absolute error (median percentage error) was
$6.2\text{\,}\mathrm{\text{Å}}$ (8.9%) for the film thickness,
$2.3\text{\,}\mathrm{\text{Å}}$ (13.3%) for the interface roughness and
$0.76\text{\times}{10}^{-6}\text{\,}{\mathrm{\text{Å}}}^{-2}$ (7.2%) for the
SLD, which is 4%, 19% and 6% higher, respectively, compared to the base model.

From this we conclude that the base model (without the added Fourier
transform) had likely already learned to implicitly extract all available
frequency information from the data and adding the Fourier components
explicitly does not lead to a better training result. Furthermore, the reason
why the results are slightly worse when the Fourier transform is added might
be attributed to the increased number of trainable parameters due to the
larger number of neurons in the model. Thus, more parameters need to be
optimized to achieve the best training result, which is generally a more
difficult task. For these reasons and the added computational requirements
during both training and inference time, we decided not to include the Fourier
transform layer into the default neural network layer of our analysis
pipeline. Nevertheless, we do not rule out that a suitable implementation of
the Fourier transform could be beneficial for certain scattering geometries.

###  VII Conclusion

We demonstrated an optimized analysis pipeline, _mlreflect_ , based on machine
learning for the automated analysis of reflectivity data. We tested our
pipeline on a large dataset of 242 XRR curves, containing in situ and ex situ
measurements of organic thin films on Si/SiOx substrates, where it showed a
performance comparable to a manually supervised least mean squares fit for
most of the data. Therefore we conclude that _mlreflect_ is a useful tool for
the automated pre-screening or even on-the-fly analysis of reflectivity data.

We also discussed that for the effective evaluation of trained machine
learning models a sufficiently large experimental dataset is necessary. Most
studies so far mainly focus on the performance of the model with regards to
simulated data and include only few, if any, experimental test data. However,
this may be misleading, since our results clearly show that the performance on
simulated data cannot easily be generalized to experimental conditions.

Furthermore, we showed the influence of possible systematic errors (such as
misalignment) on the reflectivity data and how the prediction speed of the
neural network model can be exploited to improve the overall performance by
transforming the data slightly. Our results highlight the necessity of
accounting for these differences between simulated theoretical models and real
data in order to obtain stable results.

Although the results shown here are demonstrated with systems of one layer on
a Si/SiOx substrate, the demonstrated neural network model could easily be
retrained to determine any single layer of any sample structure. While
determining multiple layers at once is in principle possible and has been
demonstrated before, this type of neural network architecture might not be
ideal to tackle this type of inverse problem with multiple solutions since
they map exactly one solution to a given input. Therefore, architectures that
yield probabilities as an output might be more suitable for multi-layer
problems.

###  VIII Acknowledgments

This research was funded by the German Federal Ministry for Science and
Education (BMBF).

We acknowledge DESY (Hamburg, Germany), a member of the Helmholtz Association
HGF, for the provision of experimental facilities. Parts of this research were
carried out at PETRA III and we would like to thank Rene Kirchhof for
assistance in using photon beamline P08. Beamtime was allocated for proposal
II-20190761. This research was also supported in part through the Maxwell
computational resources operated at DESY with the assistance of André
Rothkirch and Frank Schlünzen.

We acknowledge the Paul Scherrer Institut, Villigen, Switzerland for provision
of synchrotron radiation beamtime at the Materials Science beamline of the
SLS.

We acknowledge the European Synchrotron Radiation Facility for provision of
synchrotron radiation facilities at beamline ID10b.

We thank the Deutsche Forschungsgemeinschaft (DFG) for financial support.

Supported by the German Research Foundation through the Cluster of Excellence
“Machine Learning – New Perspectives for Science”.

### References

  * Tolan [1999] M. Tolan,  _X-ray scattering from soft-matter thin films: materials science and basic research_ , Springer Tracts in Modern Physics (Springer, Berlin, 1999).
  * Holý  _et al._ [1999] V. Holý, U. Pietsch, and T. Baumbach,  _High-Resolution X-Ray Scattering from Thin Films and Multilayers_ , Springer Tracts in Modern Physics (Springer, Berlin, 1999).
  * Braslau  _et al._ [1988] A. Braslau, P. S. Pershan, G. Swislow, B. M. Ocko, and J. Als-Nielsen, Capillary waves on the surface of simple liquids measured by x-ray reflectivity, Phys. Rev. A 38, 2457 (1988).
  * Russell [1990] T. Russell, X-ray and neutron reflectivity for the investigation of polymers, Mater Sci Rep 5, 171 (1990).
  * Neville  _et al._ [2006] F. Neville, M. Cahuzac, O. Konovalov, Y. Ishitsuka, K. Y. C. Lee, I. Kuzmenko, G. M. Kale, and D. Gidalevitz, Lipid headgroup discrimination by antimicrobial peptide ll-37: Insight into mechanism of action, Biophys. J. 90, 1275 (2006).
  * Skoda  _et al._ [2017] M. W. Skoda, B. Thomas, M. Hagreen, F. Sebastiani, and C. Pfrang, Simultaneous neutron reflectometry and infrared reflection absorption spectroscopy (IRRAS) study of mixed monolayer reactions at the air–water interface, RSC advances 7, 34208 (2017).
  * Lehmkühler  _et al._ [2008] F. Lehmkühler, M. Paulus, C. Sternemann, D. Lietz, F. Venturini, C. Gutt, and M. Tolan, The carbon dioxide-water interface at conditions of gas hydrate formation, J. Am. Chem. Soc. 131, 585 (2008).
  * Kowarik  _et al._ [2006] S. Kowarik, A. Gerlach, S. Sellner, F. Schreiber, L. Cavalcanti, and O. Konovalov, Real-time observation of structural and orientational transitions during growth of organic thin films, Phys. Rev. Lett. 96, 125504 (2006).
  * Parratt [1954] L. G. Parratt, Surface studies of solids by total reflection of x-rays, Phys. Rev. 95, 359 (1954).
  * Abelès [1950] F. Abelès, La théorie générale des couches minces, J. Phys. Radium 11, 307 (1950).
  * Heavens [1955] O. S. Heavens,  _Optical properties of thin solid films_ (London: Butterworths Scientific Publications, 1955).
  * Björck and Andersson [2007] M. Björck and G. Andersson, Genx: An extensible x-ray reflectivity refinement program utilizing differential evolution, J. Appl. Crystallogr. 40, 1174 (2007).
  * Kienzle  _et al._ [2011] P. Kienzle, J. Krycka, N. Patel, and I. Sahin, Refl1d (version 0.8.14) [computer software] (2011).
  * Nelson [2006] A. R. J. Nelson, Co-refinement of multiple-contrast neutron/x-ray reflectivity data using motofit, J. Appl. Crystallogr. 39, 273 (2006).
  * Nelson and Prescott [2019] A. R. J. Nelson and S. W. Prescott, refnx: neutron and x-ray reflectometry analysis in python, J. Appl. Crystallogr. 52, 193 (2019).
  * Danauskas  _et al._ [2008] S. M. Danauskas, D. Li, M. Meron, B. Lin, and K. Y. C. Lee, Stochastic fitting of specular x-ray reflectivity data using, J. Appl. Crystallogr. 41, 1187 (2008).
  * Gerelli [2016] Y. Gerelli, Aurore: new software for neutron reflectivity data analysis, J. Appl. Crystallogr. 49, 330 (2016).
  * Greco  _et al._ [2019] A. Greco, V. Starostin, C. Karapanagiotis, A. Hinderhofer, A. Gerlach, L. Pithan, S. Liehr, F. Schreiber, and S. Kowarik, Fast fitting of reflectivity data of growing thin films using neural networks, J. Appl. Crystallogr. 52, 1342 (2019).
  * Mironov  _et al._ [2021] D. Mironov, J. H. Durant, R. Mackenzie, and J. F. K. Cooper, Towards automated analysis for neutron reflectivity, Mach. Learn.: Sci. Technol. 2, 035006 (2021).
  * Doucet  _et al._ [2021] M. Doucet, R. K. Archibald, and W. T. Heller, Machine learning for neutron reflectometry data analysis of two-layer thin films, Mach. Learn.: Sci. Technol. 2, 035001 (2021).
  * Loaiza and Raza [2021] J. M. C. Loaiza and Z. Raza, Towards reflectivity profile inversion through artificial neural networks, Mach. Learn.: Sci. Technol. 2, 025034 (2021).
  * Greco  _et al._ [2021] A. Greco, V. Starostin, A. Hinderhofer, A. Gerlach, M. W. A. Skoda, S. Kowarik, and F. Schreiber, Neural network analysis of neutron and x-ray reflectivity data: pathological cases, performance and perspectives, Mach. Learn.: Sci. Technol. 2, 045003 (2021).
  * Andrejevic  _et al._ [2021] N. Andrejevic, Z. Chen, T. Nguyen, L. Fan, H. Heiberger, V. Lauter, L.-J. Zhou, Y.-F. Zhao, C.-Z. Chang, A. Grutter, and M. Li, Elucidating proximity magnetism through polarized neutron reflectometry and machine learning, arXiv:1410.5093 (2021), arXiv:2109.08005 [cond-mat.mtrl-sci] .
  * Abadi  _et al._ [2016] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mane, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viegas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, Tensorflow: Large-scale machine learning on heterogeneous distributed systems, arXiv:1603.04467 (2016), arXiv:1603.04467 [cs.DC] .
  * Gibaud  _et al._ [1993] A. Gibaud, G. Vignaud, and S. K. Sinha, The correction of geometrical factors in the analysis of x-ray reflectivity, Acta Cryst. A 49, 642 (1993).
  * Als-Nielsen and McMorrow [2011] J. Als-Nielsen and D. McMorrow,  _Elements of Modern X-ray Physics_, 2nd ed. (John Wiley & Sons, Ltd, Chichester, 2011).
  * Névot and Croce [1980] L. Névot and P. Croce, Characterisation of surfaces by grazing x-ray reflection, Revue de Physique Appliquée 15, 761 (1980).
  * Moré [1977] J. J. Moré, The Levenberg-Marquardt algorithm: Implementation and theory, in  _Numerical Analysis_ , Graduate Texts in Mathematics, Vol. 630, edited by G. A. Watson (Springer, New York, 1977) p. 105.
  * Smilgies  _et al._ [2005] D.-M. Smilgies, N. Boudet, B. Struth, and O. Konovalov, Troika II: a versatile beamline for the study of liquid and solid interfaces, J. Synchrotron Rad. 12, 329 (2005).
  * Seeck  _et al._ [2011] O. H. Seeck, C. Deiter, K. Pflaum, F. Bertam, A. Beerlink, H. Franz, J. Horbach, H. Schulte-Schrepping, B. M. Murphy, M. Greve, and O. Magnussen, The high-resolution diffraction beamline p08 at PETRA III, J. Synchrotron Rad. 19, 30 (2011).
  * Patterson  _et al._ [2005] B. Patterson, R. Abela, H. Auderset, Q. Chen, F. Fauth, F. Gozzo, G. Ingold, H. Kühne, M. Lange, D. Maden, D. Meister, P. Pattison, T. Schmidt, B. Schmitt, C. Schulze-Briese, M. Shi, M. Stampanoni, and P. Willmott, The materials science beamline at the swiss light source: design and realization, Nucl. Instrum. Methods Phys. Res. A 540, 42 (2005).
  * Storn and Price [1997] R. Storn and K. Price, Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces, J Global Optim 11, 341 (1997).
  * Sivia [2011] D. S. Sivia,  _Elementary scattering theory: for X-ray and neutron users_ (Oxford University Press, 2011).

## Supporting Information

###  I Optimal training noise levels

Figure 4 of the main manuscript shows the loss on the experimental dataset of
242 curves for 11 different neural network models where the training data was
modified with different amounts of uniform noise. The results show that there
seems to be an optimal noise value of about 0.3 where the loss for the
experimental data has a minimum.

An interesting question arises about how this value is related to the amount
of noise in the experimental data. To investigate this, the test data was
separated into four groups with varying amounts of noise. While the noise in
the data is not uniformly distributed, an equivalent noise level (ENL) can be
calculated by subtracting the ground truth fit from the data and taking the
absolute mean of all data points. Figure S1a shows the distribution of the ENL
across the entire dataset and how the distribution was split into the four
subsets with a different ENL. Figure S1b shows the optimal training noise (for
which the loss had a minimum) for each of the four categories as well as the
entire dataset. The error bars represent the standard deviation of five
independent training repetitions. Evidently, the ENL of the data does not seem
to have a strong influence on the optimal noise level except for the 0.4–0.5
category, where it is slightly lower. This is due to the main source of error
in the data not being statistical (e.g. Poisson noise), but rather systematic
in nature (e.g. the fit does not fully describe the data). Since the role of
the uniform noise on the training data is not to mimic the noise in the data,
but to account for these systematic deviations, the entire dataset benefits
from a similar training noise level.

However, for data with significantly higher statistical noise than our
dataset, it could be possible that optimal training noise is different.

Figure S1: (a) Distribution of the equivalent noise level (ENL) in the
experimental testing dataset of 242 XRR curves. The dataset was split into
four categories with varying ENLs to test each separately. (b) Optimal
training noise for different ENLs in the training data. The optimal level for
the entire dataset corresponds to the minimum shown in Figure 4 of the main
manuscript. The error bars represent the standard deviation of five
independent training repetitions.

###  II Detailed prediction error histograms

This section shows detailed histograms of the absolute error distribution of
each parameter with respect to the ground truth (GT) of a given parameter,
which expands on the condensed form shown in Figure 3 of the main manuscript.
Here, only the results of the full pipeline are shown (neural network + q
shift + LMS fit). Figure S2–S4 show the errors with respect to the GT
thickness, Figure S5–S7 with respect to the GT roughness and Figure S8–S10
with respect to the GT SLD.

The majority of outliers are due to ambiguous fits (e.g. featureless curves)
where multiple parameter combinations lead to a good fit. A common case are
very thin films where there are no oscillations visible in the chosen $q$
range.

Figure S2: Distribution of the absolute thickness error from the full pipeline
fit with respect to the ground truth (GT) thickness. Each dot represents a
single curve in the testing dataset. Figure S3: Distribution of the absolute
roughness error from the full pipeline fit with respect to the ground truth
(GT) thickness. Each dot represents a single curve in the testing dataset.
Figure S4: Distribution of the absolute SLD error from the full pipeline fit
with respect to the ground truth (GT) thickness. Each dot represents a single
curve in the testing dataset. Figure S5: Distribution of the absolute
thickness error from the full pipeline fit with respect to the ground truth
(GT) roughness. Each dot represents a single curve in the testing dataset.
Figure S6: Distribution of the absolute roughness error from the full pipeline
fit with respect to the ground truth (GT) roughness. Each dot represents a
single curve in the testing dataset. Figure S7: Distribution of the absolute
SLD error from the full pipeline fit with respect to the ground truth (GT)
roughness. Each dot represents a single curve in the testing dataset. Figure
S8: Distribution of the absolute thickness error from the full pipeline fit
with respect to the ground truth (GT) SLD. Each dot represents a single curve
in the testing dataset. Figure S9: Distribution of the absolute roughness
error from the full pipeline fit with respect to the ground truth (GT) SLD.
Each dot represents a single curve in the testing dataset. Figure S10:
Distribution of the absolute SLD error from the full pipeline fit with respect
to the ground truth (GT) SLD. Each dot represents a single curve in the
testing dataset.


#  Quantum Model of Bertrand Duopoly

Salman Khan, M. Ramzan , M. K. Khan
Department of Physics Quaid-i-Azam University
Islamabad 45320, Pakistan<EMAIL_ADDRESS>
###### Abstract

We present the quantum model of Bertrand duopoly and study the entanglement
behavior on the profit functions of the firms. Using the concept of optimal
response of each firm to the price of the opponent, we found only one Nash
equilibirum point for maximally entangled initial state. The very presence of
quantum entanglement in the initial state gives payoffs higher to the firms
than the classical payoffs at the Nash equilibrium. As a result the dilemma
like situation in the classical game is resolved.
PACS: 03.65.Ta; 03.65.-w; 03.67.Lx
Keywords: Quantum Bertrand duopoly; profit functions, Nash equilibria.

In economics, oligopoly refers to a market condition in which sellers are so
few that action of each seller has a measurable impact on the price and other
market factors [1]. If the number of firms competing on a commodity in the
market is just two, the oligopoly is termed as duopoly. The competitive
behavior of firms in oligopoly makes it suitable to be analyzed by using the
techniques of game theory. Cournot and Bertrand models are the two oldest and
famous oligopoly models [2, 3]. In Cournot model of oligopoly firms put
certain amount of homogeneous product simultaneously in the market and each
firm tries to maximize its payoff by assuming that the opponent firms will
keep their outputs constant. Later on Stackelberg introduced a modified form
of Cournot oligopoly in which the oligopolistic firms supply their products in
the market one after the other instead of their simultaneous moves. In
Stackelberg duopoly the firm that moves first is called leader and the other
firm is the follower [4]. In Bertrand model the oligopolistic firms compete on
price of the commodity, that is, each firm tries to maximize its payoff by
assuming that the opponent firms will not change the prices of their products.
The output and price are related by the demand curve so the firms choose one
of them to compete on leaving the other free. For a homogeneous product, if
firms choose to compete on price rather than output, the firms reach a state
of Nash equilibrium at which they charge a price equal to marginal cost. This
result is usually termed as Bertrand paradox, because practically it takes
many firms to ensure prices equal to marginal cost. One way to avoid this
situation is to allow the firms to sell differentiated products [1].

For the last one decade quantum game theorists are attempting to study
classical games in the domain of quantum mechanics [$5$-$14$]. Various quantum
protocols have been introduced in this regard and interesting results have
been obtained [15-25]. The first quantization scheme was presented by Meyer
[15] in which he quantized a simple penny flip game and showed that a quantum
player can always win against a classical player by utilizing quantum
superposition.

In this letter, we extend the classical Bertrand duopoly with differentiated
products to quantum domain by using the quantization scheme proposed by
Marinatto and Weber [17]. Our results show that the classical game becomes a
subgame of the quantum version. We found that entanglement in the initial
state of the game makes the players better off. Before presenting the
calculation of quantization scheme, we first review the classical model of the
game.

Consider two firms A and B producing their products at a constant marginal
cost $c$ such that $c<a$, where $a$ is a constant. Let $p_{1}$ and $p_{2}$ be
the prices chosen by each firm for their products, respectively. The
quantities $q_{A}$ and $q_{B}$ that each firm sells is given by the following
key assumption of the classical Bertrand duopoly model

$q_{A}=a-p_{1}+bp_{2}$ $q_{B}=a-p_{2}+bp_{1}$ (1)

where the parameter $0<b<1$ shows the amount of one firm’s product substituted
for the other firm’s product. It can be seen from Eq. (1) that more quantity
of the product is sold by the firm which has low price compare to the price
chosen by his opponent. The profit function of the two firms are given by

$u_{A}\left(p_{1},p_{2},b\right)=q_{A}\left(p_{1}-c\right)=\left(a-p_{1}+bp_{2}\right)\left(p_{1}-c\right)$
$u_{B}\left(p_{1},p_{2},b\right)=q_{B}\left(p_{2}-c\right)=\left(a-p_{2}+bp_{1}\right)\left(p_{2}-c\right)$
(2)

In Bertrand duopoly the firms are allowed to change the quantity of their
product to be put in the market and compete only in price. A firm changes the
price of its product by assuming that the opponent will keep its price
constant. Suppose that firm B has chosen $p_{2}$ as the price of his product,
the optimal response of firm A to $p_{2}$ is obtained by maximizing its profit
function with respect to its own product’s price, that is, $\partial
u_{A}/\partial p_{1}=0$, this leads to

$p_{1}=\frac{1}{2}\left(bp_{2}+a+c\right)$ (3)

Firm B response to a fixed price $p_{1}$ of firm A is obtained in a similar
way and is given by

$p_{2}=\frac{1}{2}\left(bp_{1}+a+c\right)$ (4)

Solution of Eqs.(3 and 4) lead to the following optimal price level that
defines the Nash equilibrium of the game

$p_{1}^{\ast}=p_{2}^{\ast}=\frac{a+c}{2-b}$ (5)

The profit functions of the firms at the Nash equilibrium become

$u_{A}^{\ast}=u_{B}^{\ast}=\left[\frac{a+c}{2-b}-c\right]^{2}$ (6)

From Eq. (6), we see that both firms can be made better off if they choose
higher prices, that is, the Nash equilibrium is Pareto inefficient.

To quantize the game, we consider that the game space of each firm is a two
dimensional Hilbert space of basis vector $|0\rangle$ and $|1\rangle$, that
is, the game consists of two qubits, one for each firm. The composite Hilbert
space $\mathcal{H}$ of the game is a four dimensional space which is formed as
a tensor product of the individual Hilbert spaces of the firms, that is,
$\mathcal{H=H}_{A}\mathcal{\otimes H}_{B}$, where $\mathcal{H}_{A}$ and
$\mathcal{H}_{B}$ are the Hilbert spaces of firms A and B, respectively. To
manipulate their respective qubits each firm can have only two strategies $I$
and $C$. Where $I$ is the identity operator and and $C$ is the inversion
operator also called Pauli spin flip operator. If $x$ and $1-x$ stand for the
probabilities of $I$ and $C$ that firm A applies and $y$, $1-y$, are the
probabilities that firm B applies, then the final state $\rho_{f}$ of the game
is given by [17]

$\displaystyle\rho_{f}$ $\displaystyle=$ $\displaystyle xyI_{A}\otimes I_{B}\
\rho_{i}\ I_{A}^{{\dagger}}\otimes
I_{B}^{{\dagger}}+x\left(1-y\right)I_{A}\otimes C_{B}\ \rho_{i}\
I_{A}^{{\dagger}}\otimes C_{B}^{{\dagger}}$ (7)
$\displaystyle+y\left(1-x\right)C_{A}\otimes I_{B}\ \rho_{i}\
C_{A}^{{\dagger}}\otimes I_{B}^{{\dagger}}$
$\displaystyle+\left(1-x\right)\left(1-y\right)C_{A}\otimes C_{B}\ \rho_{i}\
C_{A}^{{\dagger}}\otimes C_{B}^{{\dagger}}$

In Eq. (7) $\rho_{i}=|\psi_{i}\rangle\langle\psi_{i}|$ is the initial density
matrix with initial state $|\Psi_{i}\rangle$, which is given by

$|\psi_{i}\rangle=\cos\gamma|00\rangle+\sin\gamma|11\rangle$ (8)

where $\gamma\in\left[0,\pi\right]$ and represents the degree of entanglement
of the initial state. In Eq. (8) the first qubit corresponds to firm A and the
second qubit corresponds to firm B. The moves (prices) of the firms and the
probabilities $x$, $y$ of using the operators can be related as follows,

$x=\frac{1}{1+p_{1}},\qquad\qquad\qquad y=\frac{1}{1+p_{2}}$ (9)

where the prices $p_{1}$ and $p_{2}\in[0,\infty)$ and the probabilities $x$,
$y\in\left[0,1\right]$. By using Eqs. (7 \- 9), the nonzero elements of the
final density matrix are obtained as

$\displaystyle\rho_{11}$ $\displaystyle=$
$\displaystyle\frac{(\cos^{2}\gamma+p_{1}p_{2}\sin^{2}\gamma)}{\left(1+p_{1}\right)\left(1+p_{2}\right)}$
$\displaystyle\rho_{14}$ $\displaystyle=$
$\displaystyle\rho_{41}=\frac{(1+p_{1}p_{2})\cos\gamma\sin\gamma}{\left(1+p_{1}\right)\left(1+p_{2}\right)}$
$\displaystyle\rho_{22}$ $\displaystyle=$
$\displaystyle\frac{p_{2}\cos^{2}\gamma+p_{1}\sin^{2}\gamma}{\left(1+p_{1}\right)\left(1+p_{2}\right)}$
$\displaystyle\rho_{23}$ $\displaystyle=$
$\displaystyle\rho_{32}=\frac{(p_{1}+p_{2})\cos\gamma\sin\gamma}{\left(1+p_{1}\right)\left(1+p_{2}\right)}$
$\displaystyle\rho_{33}$ $\displaystyle=$
$\displaystyle\frac{p_{1}\cos^{2}\gamma+p_{2}\sin^{2}\gamma}{\left(1+p_{1}\right)\left(1+p_{2}\right)}$
$\displaystyle\rho_{44}$ $\displaystyle=$
$\displaystyle\frac{p_{1}p_{2}\cos^{2}\gamma+\sin^{2}\gamma}{\left(1+p_{1}\right)\left(1+p_{2}\right)}$
(10)

The payoffs of the firms can be found by the following trace operations

$\displaystyle u_{A}\left(p_{1},p_{2},b\right)$ $\displaystyle=$
$\displaystyle\mathrm{Trace}\left(U_{A}^{\mathrm{oper}}\rho_{f}\right)$
$\displaystyle u_{B}\left(p_{1},p_{2},b\right)$ $\displaystyle=$
$\displaystyle\mathrm{Trace}\left(U_{B}^{\mathrm{oper}}\rho_{f}\right)$ (11)

where $U_{A}^{\mathrm{oper}}$ and $U_{B}^{\mathrm{oper}}$ are payoffs
operators of the firms, which we define these as

$\displaystyle U_{A}^{\mathrm{oper}}$ $\displaystyle=$
$\displaystyle\frac{q_{A}}{p_{12}}\left(k_{B}\rho_{11}-\rho_{22}+\rho_{33}\right)$
$\displaystyle U_{B}^{\mathrm{oper}}$ $\displaystyle=$
$\displaystyle\frac{q_{A}}{p_{12}}\left(k_{A}\rho_{11}+\rho_{22}-\rho_{33}\right)$
(12)

where $k_{A}=p_{1}-c$, $k_{B}=p_{2}-c$ and
$p_{12}=\frac{1}{\left(1+p_{1}\right)\left(1+p_{2}\right)}$. By using Eqs. (10
\- 12), the payoffs of the firms are obtained as

$\displaystyle u_{A}\left(p_{1},p_{2},b\right)$ $\displaystyle=$
$\displaystyle(a-p_{1}+bp_{2})[k_{A}\cos^{2}\gamma+\\{p_{2}+p_{1}(-1-cp_{2}+p_{2}^{2})\\}\sin^{2}\gamma]$
$\displaystyle u_{B}\left(p_{1},p_{2},b\right)$ $\displaystyle=$
$\displaystyle(a-p_{2}+bp_{1})[k_{B}\cos^{2}\gamma+\\{p_{1}-p_{2}(1+cp_{1}-p_{1}^{2})\\}\sin^{2}\gamma]$

One can easily see from Eq. (LABEL:E13) that the classical payoffs can be
reproduced by setting $\gamma=0$ in Eq. (LABEL:E13).

We proceed similar to the classical Bertrand duopoly to find the response of
each firm to the price chosen by the opponent firm. For firm B’s price
$p_{2}$, the optimal response of firm A is obtained by maximizing its own
payoff (Eq. (LABEL:E13)) with respect to $p_{1}$. Similarly, the reaction
function of firm B to a known $p_{1}$ is obtained. These reaction functions
can be written as

$\displaystyle p_{1}$ $\displaystyle=$
$\displaystyle\frac{k_{B}[-1+p_{2}(a+bp_{2})]+[c+p_{2}+2bp_{2}-bp_{2}^{2}k_{B}+a\\{2-p_{2}k_{B}\\}]\cos
2\gamma}{(2-p_{2}k_{B})\cos 2\gamma-2p_{2}k_{B}}$ $\displaystyle p_{2}$
$\displaystyle=$
$\displaystyle\frac{k_{A}[-1+p_{1}(a+bp_{1})]+[c+p_{1}+2bp_{1}+bp_{1}^{2}k_{A}+a\\{2+p_{1}k_{A}\\}]\cos
2\gamma}{(2-p_{1}k_{A})\cos 2\gamma+2p_{1}k_{A}}$

The results of Eq. (LABEL:E14) reduce to the classical results given in Eqs.
(3 and 4) for the initially unentangled state, that leads to the classical
Nash equilibrium. This shows that the classical game is a subgame of the
quantum game.

Now, we discuss the behavior of entanglement in the initial state on the game
dynamics. It can be seen from Eq. (LABEL:E14) that the optimal responses of
the firms to a fixed price of the opponent firm, for a maximally entangled
state, are given by

$\displaystyle p_{1}$ $\displaystyle=$
$\displaystyle\frac{bp_{2}^{2}+ap_{2}-1}{2p_{2}}$ $\displaystyle p_{2}$
$\displaystyle=$ $\displaystyle\frac{bp_{1}^{2}+ap_{1}-1}{2p_{1}}$ (15)

Solving these equations, we can obtain the optimal price levels and the
corresponding payoffs of each firm. In this case the following four points are
obtained

$\displaystyle p_{1}^{\ast}(1)$ $\displaystyle=$ $\displaystyle
p_{2}^{\ast}(1)=\frac{a+\sqrt{a^{2}+4\beta}}{-2\beta}$ $\displaystyle
p_{1}^{\ast}(2)$ $\displaystyle=$ $\displaystyle
p_{2}^{\ast}(2)=\frac{2}{a+\sqrt{a^{2}+4\beta}}$ $\displaystyle
p_{1}^{\ast}(3,4)$ $\displaystyle=$
$\displaystyle\frac{2b}{a\sqrt{2+b}\left(\sqrt{2+b}\pm\gamma\right)}$
$\displaystyle p_{2}^{\ast}(3,4)$ $\displaystyle=$
$\displaystyle-\frac{1}{2b}\left[a\pm\frac{\gamma}{\sqrt{2+b}}\right]$ (16)

where the numbers in the parentheses correspond to the respective points (the
symbols $\pm$ correspond to points $3$ and $4$ respectively). To verify which
point (points) defines the Nash equilibrium of the game, we use the second
partial derivative condition. That is, for Nash equilibrium, the strategy
(point) must be the global maximum of the payoff function, that is,
$\partial^{2}u_{A(B)}/\partial p_{1(2)}^{2}<0$ and the payoff function at the
point must be higher than the payoff function at the boundary points. It can
easily be verified that this condition is satisfied only at point $1$. Hence
point $1$ defines the Nash equilibrium of the game. The payoffs of the firms
at the Nash equilibrium become

$\displaystyle u_{A}(1)$ $\displaystyle=$ $\displaystyle
u_{B}(1)=\frac{1}{4\beta^{4}}[a^{4}+2\alpha^{2}+2a^{2}b\beta+a^{3}c\beta-a\\{\left(\beta-2\right)\beta-3\\}c\beta^{2}$
(17)
$\displaystyle\qquad+\sqrt{a^{2}+4\beta}(a^{3}+2a\alpha+c\alpha^{2}+a^{2}c\beta)]$

The new parameters introduced in Eqs. (16 and 17) are defined as $\beta=b-2$,
$\alpha=2-3b+b^{2}$. The payoffs of the firms at the Nash equilibrium must be
real and positive for the entire range of substitution parameter $b$. This
condition for marginal cost $c<1.4$ is satisfied when $a\geq 3.5$. The firms’
payoffs at the other three points become

$\displaystyle u_{A}(2)$ $\displaystyle=$ $\displaystyle
u_{B}(2)=-\frac{4}{(a+\sqrt{a^{2}+4\beta})^{4}}[a^{5}c+a(-1+b)\\{(-9+5b)c-2\sqrt{a^{2}+4\beta}\\}$
$\displaystyle\qquad-a^{3}\\{(8-5b)c+\sqrt{a^{2}+4\beta}\\}+(-1+b)^{2}(-2+c\sqrt{a^{2}+4\beta})$
$\displaystyle\qquad+a^{4}(-1+c\sqrt{a^{2}+4\beta})+a^{2}\\{6-4c\sqrt{a^{2}+4\beta}+b(-4+3c\sqrt{a^{2}+4\beta})\\}]$
$\displaystyle u_{A}(3,4)$ $\displaystyle=$
$\displaystyle\frac{(1+b)^{2}(a^{2}(2+b)^{3/2}+a(2+b)(b\sqrt{2+b}c\pm\Gamma)+b(2b\sqrt{2+b}\pm
c\Gamma\left(2+b\right)))}{(2+b)^{3/2}(a(2+b)\pm\sqrt{2+b}\Gamma)^{2}}$
$\displaystyle u_{B}(3,4)$ $\displaystyle=$
$\displaystyle-\frac{(1+b)^{2}\sqrt{2+b}\left[2ac+abc-2b\pm\sqrt{2+b}\Gamma
c\right]}{4b(2+b)^{5/2}}$ (18)

where $\Gamma=\sqrt{4b^{2}+a^{2}(2+b)}$.

(a) |  |
---|---|---
Figure 1: The payoffs of the firms at the classical and quantum Nash
equilibria against the substitution parameter $b$. The values of the parameter
$a$ and the marginal cost $c$ are chosen as $3.5$ and $0.1$, respectively. The
superscripts $C$ and $Q$ of $u$ represent the classical and quantum cases,
respectively. The subscripts $A$ stands for firm $A$.

We present a quantization scheme for the Bertrand duopoly with differentiated
products. To analyze the effect of quantum entanglement on the game dynamics,
we plot the payoffs of the firms at the classical and quantum Nash equilibria
against the substitution parameter $b$ in figure ($1$). The values of
parameters $a$ and $c$ are chosen to be $3.5$ and $0.1$, respectively. The
solid line ($u_{A}^{Q}(1)$) represents quantum mechanical payoffs and the
dotted line ($u_{A}^{C}$) represents the classical payoffs of the firms. It is
clear from the figure that quantum payoffs of the firms are higher than the
classical payoffs for the entire range of substitution parameter $b$. The
maximum entanglement in the initial state of the game makes the firms better
off. In figure (2), we plot the payoffs of the firms (Eq. 18) against the
substitution parameter $b$ at the other three points which are not the Nash
equilibria.

(a) |  |
---|---|---
Figure 2: The payoffs of the firms at the second and third points as a
function of the substitution parameter $b$. The values of the parameter $a$
and the marginal cost $c$ are chosen as $3.5$ and $0.1$, respectively. The
superscript $Q$ of $u$ represent the classical and quantum cases,
respectively. The subscripts $A$ and $B$ correspond to firms $A$ and $B$
respectively. The numbers in the parentheses represent the corresponding Nash
equilibrium points.

In conclusion, we have used the Marinatto and Weber quantization scheme to
find the quantum version of Bertrand duopoly with differentiated products. We
have studied the entanglement behavior on the payoffs of the firms for a
maximally entangled initial state. We found that for large values of
substitution parameter $b,$ both firms can achieve significantly higher
payoffs as compared to the classical payoffs. Furthermore, for maximally
entangled state the quantum payoffs are higher than the classical payoffs for
the entire range of substitution parameter and is the best situation for both
firms. Thus, the dilemma-like situation in the classical Bertrand duopoly game
is resolved.

Acknowledgment

One of the authors (Salman Khan) is thankful to World Federation of Scientists
for partially supporting this work under the National Scholarship Program for
Pakistan

## References

  * [1] Gibbons R Game Theory for applied Economists Princeton Univ. Press Princeton, NJ, 1992: Bierman H.S, Fernandez L Game Theory with Economic Applications, 2nd Edition, Addison - Wesley, Reading MA 1998
  * [2] Cournot A 1897 Researches into the Mathematical Principles of the Theory of Wealth, Macmillan Co., New York
  * [3] Bertrand J, Savants J 1883 67 499
  * [4] Stackelberg H von 1934 Marktform und Gleichgewicht (Julius Springer Vienna)
  * [5] Flitney A.P, Ng J, Abbott D 2002 Physica A 314 35, doi:10.1016/S0378-4371(02)01084-1
  * [6] Ramzan M et al 2008 J. Phys. A Math. Theor. 41 055307, doi:10.1088/1751-8113/41/5/055307
  * [7] D’ Ariano G M 2002 Quantum Information and Computation 2 355
  * [8] Flitney A.P, Abbott D 2002 Phys. Rev. 65 062318, doi:10.1103/PhysRevA.65.062318
  * [9] Iqbal A, Cheo T, Abbott D 2008 Phys. Lett. A 372 6564, doi:10.1016/j.physleta.2008.09.026
  * [10] Zhu X, Kaung L M 2008 Commun. Theor. Phys (china) 49 111
  * [11] Ramzan M and Khan M K 2008 J. Phys. A: Math. Theor. 41 435302, doi:10.1088/1751-8113/41/43/435302
  * [12] Ramzan M, Khan M.K 2009 J. Phys. A, Math. Theor. 42 025301, doi:10.1088/1751-8113/42/2/025301
  * [13] Zhu X, Kaung L M 2007 J. Phys. A 40 7729, doi:10.1088/1751-8113/40/27/021
  * [14] Salman Khan et al 2010 Int. J. Theor. Phys.49 31, doi:10.1007/s10773-009-0175-y
  * [15] Meyer D A 1999 Phys. Rev. Lett. 82 1052, doi:10.1103/PhysRevLett.82.1052
  * [16] Eisert J et al 1999 Phys. Rev. Lett. 83 3077, doi:10.1103/PhysRevLett.83.3077
  * [17] Marinatto L and Weber T 2000 Phys. Lett. A 272 291, doi:10.1016/S0375-9601(00)00441-2
  * [18] Li H, Du J, Massar S 2002 Phys. Lett. A 306 73, doi:10.1016/S0375-9601(02)01628-6
  * [19] Lo C.F, Kiang D 2004 Phys. Lett. A 321 94, doi:10.1016/j.physleta.2003.12.013
  * [20] Iqbal A and Abbott D 2009 arXiv:quant-ph/0909.3369
  * [21] Iqbal A and Toor A H 2002 Phys. Rev. A 65 052328, doi:10.1103/PhysRevA.65.052328
  * [22] Lo C.F, Kiang D 2003 Phys. Lett. 318 333, doi:10.1016/j.physleta.2003.09.047
  * [23] Benjamin S.C, Hayden P.M 2001 Phy. Rev. A 64 030301, doi:10.1103/PhysRevA.64.030301
  * [24] Lo C.F, Kiang D 2005 Phys. Lett. A 346 65, doi:10.1016/j.physleta.2005.07.055
  * [25] Gan Qin et al 2005 J. Phys. A: Math. Gen. 38 4247, doi:10.1088/0305-4470/38/19/013