# Angle-resolved optically detected magnetic resonance as a tool for strain
determination in nanostructures

A. Bogucki<EMAIL_ADDRESS>Institute of Experimental Physics,
Faculty of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland
M. Goryca  Institute of Experimental Physics, Faculty of Physics, University
of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  A. Łopion  Institute of
Experimental Physics, Faculty of Physics, University of Warsaw, Pasteura 5,
02-093 Warsaw, Poland  W. Pacuski  Institute of Experimental Physics, Faculty
of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  K.
Połczyńska  Institute of Experimental Physics, Faculty of Physics, University
of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  J. Domagała  Institute of
Physics, Polish Academy of Sciences, Aleja Lotników 32/36, 02-668 Warsaw,
Poland  M. Tokarczyk  Institute of Experimental Physics, Faculty of Physics,
University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  T. Fąs  Institute of
Experimental Physics, Faculty of Physics, University of Warsaw, Pasteura 5,
02-093 Warsaw, Poland  A. Golnik  Institute of Experimental Physics, Faculty
of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  P.
Kossacki  Institute of Experimental Physics, Faculty of Physics, University of
Warsaw, Pasteura 5, 02-093 Warsaw, Poland

###### Abstract

In this paper, we apply the angle-resolved Optically Detected Magnetic
Resonance (ODMR) technique to study series of strained (Cd, Mn)Te/(Cd, Mg)Te
quantum wells (QWs) produced by molecular beam epitaxy. By analyzing
characteristic features of ODMR angular scans, we determine strain-induced
axial-symmetry spin Hamiltonian parameter $D$ with neV precision. Furthermore,
we use low-temperature optical reflectivity measurements and X-ray diffraction
scans to evaluate the local strain present in QW material. In our analysis, we
take into account different thermal expansion coefficients of GaAs substrate
and CdTe buffer. The additional deformation due to the thermal expansion
effects has the same magnitude as deformation origination from the different
compositions of the samples. Based on the evaluated deformations and values of
strain-induced axial-symmetry spin Hamiltonian parameter $D$, we find strain
spin-lattice coefficient ${G_{11}=(72.2\pm 1.9)}$ neV for Mn2+ in CdTe and
shear deformation potential ${b=(-0.94\pm 0.11)}$ eV for CdTe.

###### pacs:

76.70.Hb,78.67.De,78.55.Et,78.66.Hf,78.70.Gq ,75.50.Pp,71.70.Fk,72.25.Rb

##  I Introduction

One of the critical factors influencing the electronic and optoelectronic
devices’ performance is local strain distribution. For example, the spin
relaxation time – a crucial parameter for potential spintronic devices –
strongly depends on spin-lattice coupling and local strain distribution. The
latter is particularly non-trivial in complex nanostructures composed of many
different materials, e.g., in quantum wells (QW). Such structures can be
produced with several different growth techniques like metalorganic vapor-
phase epitaxy (MOVPE) or molecular beam epitaxy (MBE). Historically, the
development of semiconductors’ growth methods was accompanied by independent
improvement of materials characterization techniques, which resulted in a
comprehensive knowledge foundation for future research and advanced
applications. Naturally, new characterization techniques combined with new
growth methods offer fresh insights into semiconductor physics and often shed
new light on some previously determined material parameters. For example, this
was the case of GaAs crystals widely used as substrates for complex
semiconductor structures. In 1959 Kolm et al. Kolm  _et al._ (1957) used X-ray
diffraction to determine the lattice constant of the GaAs crystals produced
with the Bridgman–Stockbarger method. However, crystals grown nowadays with
the use of other techniques like Horizontal Bridgman (HB), Liquid Encapsulated
Czochralski (LEC), Vertical Gradient Freeze (VGF), or epitaxial methods
resulting in different concentrations of impurities systematically exhibit
different lattice constant Bassignana  _et al._ (1997); Usuda and Ando (1996).

Finding a particular material parameter for materials obtained by a specific
growth method is critical for real multilayer devices. Determining some
structural parameters is particularly challenging in small-volume structures
(e.g. single QWs). A small amount of material often excludes many
characterization techniques like X-ray diffraction which is a standard method
of measuring the strain in semiconductor nanostructures. A possible solution
to this issue is an incorporation of a small amount (below 1%) of magnetic
ions into the structure and determination of local strain exploiting the
coupling between those ions and the strain of the crystal lattice. In many
cases (like multiple QWs structures), determining the spin Hamiltonian
parameters describing such coupling may be done with electron paramagnetic
resonance (EPR) techniques Qazzaz  _et al._ (1995). However, in thin
structures, the number of spins is insufficient to produce a detectable EPR
signal. This problem can be overcome by using the optically detected magnetic
resonance (ODMR) technique which exploits strong exchange interaction between
magnetic ions and photo-generated carriers and facilitates detection of the
paramagnetic resonance of significantly smaller amount of magnetic ions.

Figure 1: The schematic overview of this work. There are three main branches
related to the experimental techniques used in this paper: angular resolved
ODMR measurements, reflectivity measurements, and XRD measurements. By
combining the results of these methods we determine the spin-lattice coupling
parameter $G_{11}$ and deformation potential parameter $b$ for Mn2+ ion in
CdTe.

In this paper, we present a comprehensive analysis of the local strain present
in series of (Cd, Mn)Te/(Cd, Mg)Te QWs. We combine data obtained from angle-
resolved ODMR, low-temperature reflectance measurements, and X-ray diffraction
scans. The scheme presented in figure 1 shows how we address problems
formulated in the previous paragraph. We design a series of (Cd, Mn)Te/(Cd,
Mg)Te quantum wells samples fabricated by MBE method that have different
deformation in the QW layer by design. The QW material is doped with a small
amount of manganese ($<1$%) which exhibits paramagnetic resonance. Using
angle-resolved ODMR, we determine strain-induced axial-symmetry spin
Hamiltonian parameter ($D$) with high accuracy. Independently, we perform XRD
measurements on the least and the most strained samples – as only for these
samples the XRD method yields reliable deformation value. After considering
additional deformation that originates from different thermal expansion
coefficients of the substrate and QW material, we deduce deformation values
for all samples. To check if other samples also follow the design scheme, we
perform low-temperature reflectivity measurements. We determine energy
splitting between light-hole (LH) exciton transition and heavy-hole (HH)
exciton transition from reflectivity spectra. We find that the HH-LH splitting
for the whole series of samples exhibits monotonic behavior - thus validating
the designed distribution of deformation. As the ODMR measurements are taken
at cryogenic temperatures, we evaluate CdTe elastic stiffness constants
$C_{11}$ and $C_{12}$ at 1.6 K from available literature data. Then from
obtained deformation values and corresponding strain-induced axial-symmetry
spin Hamiltonian parameter $D$ combined with elastic constants, we calculate
the spin-lattice coupling coefficient $G_{11}$ for manganese in cadmium
telluride. The spin-lattice coupling coefficient $G_{11}$ value obtained in
this work diverges from previously reported in ref. Causa and Tovar, 1980.
Furthermore, we use HH-LH splittings and deformation data to find deformation
potential $b$. The obtained value of $b$ is similar to ref’s Merle d’Aubigné
_et al._ , 1990 one but different from values reported in refs. Thomas, 1961;
Mathieu  _et al._ , 1988; Peyla  _et al._ , 1992.

Our findings could be especially useful in all studies involving spin-related
phenomena in CdTe-based MBE-grown nanostructures and devices including
multiple-QWs and quantum dots Cherbunin  _et al._ (2020); Lafuente-Sampietro
_et al._ (2017); Dinu  _et al._ (2017); Moldoveanu  _et al._ (2016); Lafuente-
Sampietro  _et al._ (2016, 2015); Goryca  _et al._ (2015); Cronenberger  _et
al._ (2015); Kobak  _et al._ (2014); Varghese  _et al._ (2014); Besombes and
Boukari (2014); Goryca  _et al._ (2014); Cronenberger  _et al._ (2013); Jamet
_et al._ (2013); Chen  _et al._ (2013); Besombes  _et al._ (2012); Goryca  _et
al._ (2009a, b); Le Gall  _et al._ (2009).

##  II Experiment and results

###  II.1 Samples

The samples containing QWs used in this work were produced by the Molecular
Beam Epitaxy (MBE) technique. The scheme of samples is presented in figure 2a.
In order to obtain structures with a different strain, we have changed the
content of magnesium in the buffer layer and the barrier layer (see table 1).
The magnesium content was deduced based on magnesium fluxes and calibration
obtained by the reflectance measurements Hartmann  _et al._ (1996); LeBlanc
_et al._ (2017). The nominal width of (Cd, Mn)Te QW was 10 nm. The value of
the QW width was chosen to ensure sufficient confinement together with narrow
excitonic features in the reflectivity spectra, as well as it is well below
the critical thickness of the CdTe lattice relaxation Cibert  _et al._ (1990,
1991). The manganese content was determined by magnetooptical measurements of
the giant Zeeman splitting and fitting the modified Brillouin functionGaj  _et
al._ (1994). The thickness of the (Cd, Mg)Te barriers was 50 nm, the (Cd,
Mg)Te buffer layer was 2000 nm (far above the critical thicknesses of the (Cd,
Mg)TeWaag  _et al._ (1993)) and the additional CdTe layer below the buffer was
4000 nm thick. This additional layer ensured the isolation from the GaAs
substrate. To confirm that we obtained designed strain in produced samples we
have performed reflectivity measurements (see figure 2b). The identification
of the spectrum features was based on standard magneto-optical measurements
(not shown)Kossacki  _et al._ (1997); Kossacki (2003); Kossacki  _et al._
(2004). The energy splitting between light-hole exciton and heavy-hole exciton
visible in spectroscopic measurements is a linear function of the strain
present in the QW layer Merle d’Aubigné  _et al._ (1990); Allègre  _et al._
(1990); Boley  _et al._ (1993); Bhattacharjee and Fishman (1998). Therefore,
as the samples are similarly designed, the light-hole heavy-hole splitting
gives a hint about the existing strain.

Figure 2: (Color online) (a) An schematic of a representative sample structure used in the whole sample series used in this article. The buffer layer was made of (Cd, Mg)Te where magnesium content was varying from 0% up to approx. 30%. Barriers were made of (Cd, Mg)Te while quantum well (QW) mas made of (Cd, Mn)Te with manganese concentration below 0.5%. (b) Reflectivity spectra of two example samples. Sample UW1030 has a higher concentration of magnesium in the buffer layer than UW1031 and thus exhibits larger strain which results in larger energy splitting between heavy-hole exciton (HH) and light-hole exciton (LH). Table 1: Composition of samples used in this article, deformation present in QW calculated for 300 K without correction for different thermal expansion coefficients of CdTe and GaAs, deformation present in QW calculated for 1.6 K present in QW layer with thermal-expansion-coefficients-difference correction. | buffer | barrier | QW | ${\varepsilon^{T=300~{}K}_{\parallel\mathrm{~{}w/o~{}corr.}}}$ |  ${\varepsilon^{T=1.6~{}K}_{\parallel\mathrm{~{}corr.}}}$
---|---|---|---|---|---
sample no. | Mg (%) | Mg (%) | Mn (%) | (‰) | (‰)
UW1029 | $30.7$ | $30.7$ | $0.30$ | $-2.870$ | $-3.512$
UW1030 | $21.2$ | $21.2$ | $0.31$ | $-1.959$ | $-2.601$
UW1028 | $16.4$ | $16.4$ | $0.14$ | $-1.538$ | $-2.180$
UW0677 | $15.2$ | $15.2$ | $0.26$ | $-1.396$ | $-2.038$
UW1050 | $8.5$ | $16.4$ | $0.15$ | $-0.780$ | $-1.422$
UW1031 | $0$ | $21.2$ | $0.26$ | $0.058$ | $-0.584$
UW0676 | $0$ | $17.3$ | $0.27$ | $0.060$ | $-0.582$
Figure 3: The energy difference between the lowest light-hole excitonic state
and lowest heavy-hole excitonic state ($\Delta$LHHH splitting) plotted versus
magnesium content in the buffer layer for all samples from table 1 confirms
that the growth process was fully controlled.

Figure 3 shows heavy-hole exciton light-hole exciton splitting obtained as a
function of magnesium content in the buffer layer for all samples. The
monotonic behavior of extracted data confirms that samples follow the design –
the higher magnesium content, the larger light-hole heavy-hole splitting.

###  II.2 XRD measurements

In order to independently calibrate deformation, we have performed room-
temperature X-ray diffraction (XRD) measurements of the two extreme samples:
UW1029 (highest strain) and UW1031 (lowest strain) – see figure 4. By
analyzing the 004 and the ${}\mkern 1.0mu\overline{\mkern-1.0mu3}{}\mkern
1.0mu\overline{\mkern-1.0mu3}5$ reflections, we extract lattice constants of
CdTe separation layer (about a$=6.481$ Å) and the strained layers. Throughout
this work, for the calculations we use the following lattice constants of pure
materials at room temperature: a${}_{\mathrm{CdTe}}=6.481$ Å,
a${}_{\mathrm{MnTe}}=6.337$ Å, a${}_{\mathrm{MgTe}}=6.419$ Å and
a${}_{\mathrm{GaAs}}=5.6535$ Å. Obtained results agree with the sample design.

Figure 4: (Color online) Upper panel shows XRD maps with 004 (upper left) and
${}\mkern 1.0mu\overline{\mkern-1.0mu3}{}\mkern
1.0mu\overline{\mkern-1.0mu3}5$ (upper right) reflections obtained for high-
strain sample (UW1029). There are two strong peaks attributed to CdTe and to
(Cd, Mg)Te buffer layer. Analysis of ${}\mkern
1.0mu\overline{\mkern-1.0mu3}{}\mkern 1.0mu\overline{\mkern-1.0mu3}5$ map
confirms that strain existing in the QW is induced by buffer layer, as weak
satellite peaks are present along vertical line meaning pseudomorphic growth
with buffer lattice constant. The lower panel presents maps from the same
region of reciprocal space (lower left panel the 004 reflection, lower right
panel for ${}\mkern 1.0mu\overline{\mkern-1.0mu3}{}\mkern
1.0mu\overline{\mkern-1.0mu3}5$ reflection) for the low-strain sample
(UW1031). All peaks present in ${}\mkern 1.0mu\overline{\mkern-1.0mu3}{}\mkern
1.0mu\overline{\mkern-1.0mu3}5$-map appear along the horizontal line above
CdTe peak, confirming the pseudomorphic character of growth with CdTe lattice
constant.

In XRD maps of the high-strain sample (UW1029), there are two strong
reflections – one corresponding to 4 $\upmu$m CdTe layer, and the second one
to 2 $\upmu$m (Cd, Mg)Te buffer layer. Positions of the reflections give the
lattice constants of CdTe and buffer at room temperature (where ∥ and ⟂ denote
directions parallel and perpendicular to the samples surface, respectively):
$a^{\textrm{CdTe}}_{\parallel}=(6.4807\pm 0.0006)$ Å,
$a^{\textrm{CdTe}}_{\perp}=(6.4862\pm 0.0001)$ Å,
$a^{\textrm{buff}}_{\parallel}=(6.4616\pm 0.0003)$ Å,
$a^{\textrm{buff}}_{\perp}=(6.4595\pm 0.0001)$ Å. From this values and using a
formula

$a_{\textrm{relax}}=(C_{11}a_{\perp}+2C_{12}a_{\parallel})/(C_{11}+2C_{12})$

we can calculate lattice constant of relaxed material and finally obtain
deformation of the buffer layer
$\varepsilon_{\parallel}=-0.3565\text{\textperthousand}$. The fact that
nominally relaxed material has nonzero deformation is explained in following
sections. The $C_{12}$ and $C_{11}$ are elastic stiffness constants of CdTe.
More details about evaluation of elastic stifness constants at different
temperatures are provided in the Appendix. The analysis of ${}\mkern
1.0mu\overline{\mkern-1.0mu3}{}\mkern 1.0mu\overline{\mkern-1.0mu3}5$
reflection shows that the buffer layer is not fully relaxed. The relaxation
degree of the buffer layer in UW1029 sample is $95.6\%$.

The situation is easier to interpret in the case of the low-strain sample
(UW1031) where on top of the separation layer made of CdTe there is a buffer
layer also made from CdTe. In principle, the lattice constant of the thin,
strained barrier (thickness well below critical thickness) should be matched
to the lattice constant of the thick CdTe buffer. Therefore lattice constant
of the QW should be close to the lattice constant of CdTe. The obtained values
of lattice constants for the UW1031 sample from XRD measurements are:
$a^{\textrm{CdTe+buffer}}_{\parallel}=(6.4807\pm 0.0006)$ Å,
$a^{\textrm{CdTe+buffer}}_{\perp}=(6.4867\pm 0.0001)$ Å,
$a^{\textrm{2nd}}_{\parallel}=(6.4797\pm 0.0008)$ Å,
$a^{\textrm{2nd}}_{\perp}=(6.4549\pm 0.0001)$ Å. From this, we can calculate
$\varepsilon_{\parallel}=-0.389\text{\textperthousand}$ which in the case of
the UW1031 sample is also the deformation present in QW at room temperature.

The influence of temperature is critical in the correct determination of
deformation in heterogenic structures Horning and Staudenmann (1987); Tatsuoka
_et al._ (1989, 1990, 1991) which is the case of our sample, as the used
materials exhibit different temperature expansion coefficients. There are
three distinctive temperature points that should be considered in the context
of this work. i) The growth temperature at which CdTe layers were deposited in
MBE. We have estimated this temperature to be 552 K – see Appendix. ii) The
room temperature ($300$ K) at which XRD measurements were taken. During the
cooling after the growth process, an additional compressive deformation builds
up in the QW layer, as the GaAs substrate shrinks more than CdTe during the
cooling. So at room temperature there is an additional deformation originating
from the different shrinking of the QW layer and the substrate layer that is
visible in XRD measurements:
$\varepsilon^{\textrm{552K}\rightarrow\textrm{300K}}_{\parallel}=-0.389\text{\textperthousand}$.
iii) The pumped liquid helium temperature ($\approx 1.6$ K) at which
reflectivity and ODMR measurements were performed. The GaAs temperature
expansion coefficient differs from that of the CdTe (see figure 13 in
Appendix). Moreover, the magnitude of that difference changes as the
temperature decreases finally resulting in an additional compression which
gives contribution to the deformation
$\varepsilon^{\textrm{300K}\rightarrow\textrm{1.6K}}_{\parallel}=-0.253\text{\textperthousand}$.
Taking into account all of those contributions up together with starting value
$\varepsilon^{T=300~{}K}_{\parallel\mathrm{~{}w/o~{}corr.}}=0.058\text{\textperthousand}$
we finally obtain deformation in the CdTe QW layer at 1.6 K for the UW1031
sample $\varepsilon_{\parallel}=-0.584\text{\textperthousand}$.

Now, using above calibration and taking into account an additional
temperature-originating deformation it is possible to present light-hole-
heavy-hole splitting as a function of deformation (figure 6). We obtain band
shear deformation potential ${b=(-0.94\pm 0.11)}$ eV using the formula Chuang
(2009):

$\Delta\mathrm{LHHH}=2b\left(1+\frac{2C_{12}}{C_{11}}\right)\varepsilon_{\parallel}+\mathrm{const.}$
(1)

The shear deformation potential obtained in this work is close to the value
$(-1.05\pm 0.01)$ eV obtained by Merle d‘Aubigné et al. Merle d’Aubigné  _et
al._ (1990). The other experimentally determined value is $b=1.24$ eV by
Thomas and coworkers Thomas (1961). Other values available are given by Peyla
et al. Peyla  _et al._ (1992), and by Mathieu et al. $-1.4$ eV Mathieu  _et
al._ (1988). However, the two latter ones are based on refs. Thomas (1961);
Merle d’Aubigné  _et al._ (1990). Namely, the value given by Peyla et al.
Peyla  _et al._ (1992) is the mean of values given by Thomas and Merle. The
results of calculated light-hole heavy-hole splitting obtained from numerical
solving the finite QW model are marked in Figure 6 with black crosses. More
details about QWs simulations are in Appendix.

Figure 5: (Color online) The $2\theta$ XRD scan of the 004 reflection for low-
strain (UW1031, black line) and high-strain (UW1029, blue line) samples.
Lattice constants calculated from this measurement confirm that obtained
strain agrees with the designed one.  Figure 6: The energy difference between
the lowest light-hole excitonic state and lowest heavy-hole excitonic state
($\Delta$LHHH splitting) plotted as a function of deformation
$\varepsilon_{\parallel}$ at $1.6$ K for analyzed quantum well samples. The
empty circles represent measured data, the black crosses are results of
numerical simulations described in Appendix II.2. The black line is a result
of a linear fit to the data points. The shear deformation potential
${b=(-0.94\pm 0.11)}$ eV is calculated from the slope of an obtained line
(equation 1).

###  II.3 Optically Detected Magnetic Resonance measurements

Optically Detected Magnetic Resonance is a technique that exploits the fact
that the optical properties of the studied material change when paramagnetic
resonance occurs. The microwave absorption at electron paramagnetic resonance
(EPR) frequency leads to an increase of Mn-system temperature which results in
a decrease of optically-detected magnetization van Gisbergen  _et al._ (1993);
Ivanov  _et al._ (2001); Godlewski  _et al._ (2002); Khachapuridze  _et al._
(2002); Godlewski  _et al._ (2004); Smith  _et al._ (2005); Godlewski  _et
al._ (2008); Baranov  _et al._ (2009); Debus  _et al._ (2016); O. Tolmachev
_et al._ (2020); Shornikova  _et al._ (2020). The change in the optical
properties of the sample can be detected with various experimental techniques.
For example, previous ODMR studies of Diluted Magnetic Semiconductors were
based on changes in Faraday rotation Komarov  _et al._ (1977), in the
amplitude of photoluminescence or in the spectral shift of photoluminescence
line Ivanov  _et al._ (2008); Ivanov and Godlewski (2010); Sadowski  _et al._
(2003); Byszewski  _et al._ (2004); Tolmachev  _et al._ (2010, 2012); Gurin
_et al._ (2015). In this work, we use the energy position of the neutral
exciton emission line (X) from photoexcited QW to detect local magnetization
in the QW layer (fig. 7). For optical excitation, we use $\lambda=647$ nm
laser (with a spot diameter of approx. $100\>\upmu$m) in a standard
photoluminescence experimental setup. The sample is placed in the optical
cryostat equipped with two perpendicular pairs of split superconducting coils.
Thus it is possible to obtain up to $3$ T of a magnetic field in any direction
in a given plane. In particular, it is possible to continuously change the
magnetic field from Faraday to Voigt configuration (scan over the out-of-plane
angle $\theta$ in figure 9) or to rotate the magnetic field within the sample
plane (scan over the in-plane angle $\varphi$ in figure 9).

Figure 7: (Color online) (a) Photoluminescence (PL) spectra of (Cd, Mn)Te/(Cd,
Mg)Te quantum well (QW) collected at $1.6$ K at a magnetic field of $0.545$ T
with microwave radiation ($15.3$ GHz) turned on (red curve) and off (black
curve). (b) Map of the normalized photoluminescence spectra of the same QW
measured for various magnetic fields with present microwave radiation. The
energetic position of the PL spectrum follows the modified Brillouin function,
except the points at $0.545$ T, where the paramagnetic resonance of Mn2+ ions
placed in the QW occurs. As a result of the resonance, the giant Zeeman
splitting rapidly decreases in those points. Figure 8: Optically Detected
Magnetic Resonance (ODMR) signal of (Cd, Mn)Te/(Zn, Mg)Te quantum well
obtained as a difference between the energetic position of the quantum well
excitonic line measured with microwave radiation and expected position of the
line calculated from fitted modified Brillouin function without microwave
radiation. On the the horizontal axis, we used the relative magnetic field – a
magnetic field with a subtracted value of the resonance magnetic field for a
given frequency of microwave radiation. The multiple features visible on the
spectrum are originating from different terms of manganese Hamiltonian
hyperfine structure interaction. Figure 9: (Color online) Angles definition of
the magnetic field directions. The sample edge is along $\langle$110$\rangle$
direction. The out-of-plane angle $\theta$ is the angle between the samples’
growth axis and samples’ surface. For example, $\theta=0^{\circ}$ means
magnetic field along growth axis in Faraday configuration, whereas
$\theta=90^{\circ}$ means in-plane magnetic field. The in-plane angle
$\varphi$ is the angle between $\langle$100$\rangle$ direction and direction
of the magnetic field. The angle $\varphi=45^{\circ}$ means that the in-plane
magnetic field was along $\langle$110$\rangle$ direction.

The representative photoluminescence spectrum (PL) for samples used in this
work is presented in figure 7a. The PL spectrum consists of two emission lines
— at lower energies there is a line related to the charged exciton (CX) and at
higher energies, to the neutral exciton (X)Kossacki  _et al._ (2004). During
the scan in the magnetic field in Faraday configuration (magnetic field
perpendicular to the surface of the sample and parallel to the optical axis)
energetic position of X line can be described with the modified Brillouin
function Gaj  _et al._ (1994) (giant Zeeman effect). However, in the presence
of microwave radiation, a paramagnetic resonance of manganese ions occurs, at
the resonant magnetic field, which results in a rapid decrease of giant Zeeman
splitting (see figure 7b with its caption). Comparison of PL spectra measured
for a resonance magnetic field with and without microwave radiation are
presented in figure 7a. The difference in the energy position of X line for
these two cases is called further in the text ”the ODMR amplitude/signal“. It
is important to note that too strong laser excitation can alter the optical
properties of QW resulting in a change of the position of the X line. To avoid
this effect and ensure that the laser excitation does not influence an
amplitude of the Giant Zeeman splitting Keller  _et al._ (2001); Seufert  _et
al._ (2001); Yakovlev  _et al._ (2004); Hundt  _et al._ (2005) we use the
following procedure to determine an optimal laser excitation power. The
magnetic field is set to the value of paramagnetic resonance with microwave
radiation turned off. The series of PL spectra are measured for decreasing
laser power. The optimal laser power is a power that does not cause an energy
shift of the neutral exciton emission line. In this work, the optimal laser
power corresponds to the excitation density power $\rho<0.005$ W/cm2, which is
even smaller than the low-power regime value presented in ref. Keller  _et
al._ , 2001.

An example of representative ODMR signal as a function of magnetic field for
fixed microwave frequency is presented in figure 8. The detailed features
(multiple lines) visible in the ODMR spectrum originate from interactions of
manganese ion Mn2+ and CdTe material in QW and can be described with spin
Hamiltonian Qazzaz  _et al._ (1995); Abragam and Bleaney (2012):

$\displaystyle\widehat{H}$ $\displaystyle=$ $\displaystyle
g_{\mbox{\tiny{Mn}}}\mu_{B}\mathbf{\hat{B}\hat{S}}+A\mathbf{\hat{I}\hat{S}}+D\left[\hat{S}_{z}^{2}-\frac{S(S+1)}{3}\right]+{}$
$\displaystyle{}+E\;R_{z}(\alpha)\left(\hat{S}_{x}^{2}-\hat{S}_{y}^{2}\right)R_{z}^{-1}(\alpha)+$
$\displaystyle{}+\frac{a}{6}\left[\hat{S}_{x}^{4}+\hat{S}_{y}^{4}+\hat{S}_{z}^{4}-\frac{S(S+1)(3S^{2}+3S-1)}{5}\right].$

The first term of the Hamiltonian is responsible for Zeeman splitting (where
$g_{\mbox{\tiny{Mn}}}$ is g-factor of manganese ion). The second one is a
hyperfine coupling between the electronic spin and nuclear spin and it results
in splitting the ODMR spectrum into $6$ lines as manganese $2+$ ion has
electronic spin $S=5/2$ and nuclear spin $I=5/2$. The third and fourth term
comes from electronic quadrupole fine structure present in less than cubic
symmetry (strained QW) and can be written in a general form as
$\hat{S}\mathbf{\hat{D}}\hat{S}$ where $\mathbf{\hat{D}}$ is a tensor.
However, it can be reduced to only two parameters $D$ and $E$ – the axial and
rhombic zero-field splitting (ZFS) parameters ($D$ is also called the strain-
induced axial-symmetry parameter)Abragam and Bleaney (2012). The
$R_{z}(\alpha)$ allows for the rotation of the reference frame around the
growth axis of the QW. The last term describes the crystal-field splitting
where $a$ is the zero-field fine structure splitting parameter for the
unstrained CdTe. As a consequence of the above Hamiltonian, the microwave
absorption spectrum of Mn2+ ion consists of $30$ transitions ($5$ electronic
spin transitions, each split into the sextet due to hyperfine interaction).

###  II.4 Angular resolved ODMR

Angular resolved ODMR measurements can give precise information about spin
Hamiltonian parameters $D$ and $E$. While single ODMR measurement in Faraday
configuration should in principle give information about the $D$ parameter,
the angular ODMR can provide much more detailed information and increase
experimental accuracy. The left panels of figure 11 and 11 show angular ODMR
measurements of example high-strain and low-strain samples respectively.

Figure 10: (Color online) The left panel shows measured angular maps of the
ODMR signal for the exemplary high-strain sample (UW0677). The intensity of
the ODMR signal is coded in the brightness of the map. Upper left part (in-
plane scan) was obtained for fixed out-of-plane angle $\theta=90^{\circ}$. A
lower left part (out-of-plane) scan was obtained for fixed in-plane angle
$\varphi=45^{\circ}$. The right part presents results of simulations obtained
by numerical solving of the spin Hamiltonian presented in equation II.3.
Fitted spin Hamiltonian parameters are $D=(497\pm 3)$ neV and $E=(-8\pm 13)$
neV.

Figure 11: (Color online) The left panel shows measured angular maps of the
ODMR signal for the exemplary low-strain sample (UW0676). The intensity of the
ODMR signal is coded in the brightness of the map. The upper left part (in-
plane scan) was obtained for fixed out-of-plane angle $\theta=90^{\circ}$. A
lower left part (out-of-plane) scan was obtained for fixed in-plane angle
$\varphi=45^{\circ}$. The right part presents results of simulations obtained
by numerical solving of the spin Hamiltonian presented in equation II.3.
Fitted spin Hamiltonian parameters are $D=(155\pm 3)$ neV and $E=(0\pm 25)$
neV.

In the lower part of the plot, the angular scan from the Faraday configuration
(magnetic field direction $\theta=0^{\circ}$) to the Voigt configuration
(magnetic field direction $\theta=90^{\circ}$) is presented. The upper part of
the plot shows an in-plane ODMR angular scan starting from
$\langle$110$\rangle$ direction (parallel to the samples’ edge) for
$\varphi=45^{\circ}$. The whole ODMR angular map resents a set of
characteristic features in which location on the map allows for the
determination of spin Hamiltonian parameters by comparison with similar
numerically simulated maps (see right panels).

The first characteristic feature is the width of the ODMR spectrum in Faraday
configuration. In this configuration, a transition that is visible at the
lowest magnetic field corresponds to the change of the electron spin
projection from $S_{z}=+3/2$ to $S_{z}=+5/2$ with nuclear spin projection
$I_{z}=+5/2$. Similarly, the transition at the highest magnetic field
corresponds to the change of the electron spin projection from $S_{z}=-5/2$ to
$S_{z}=-3/2$ with nuclear spin projection $I_{z}=-5/2$. The measure of
splitting between the outermost lines in units of the magnetic field directly
corresponds to the value of the $D$ parameter in the spin Hamiltonian. At high
magnetic field approximation – assuming that eigenstates are the same as the
eigenstates of the $z$ components of the spin operators, the spin Hamiltonian
can be analytically solved. In such a case the $D$ parameter can be calculated
as

$D=(g_{\mbox{\tiny{Mn}}}\mu_{B}\Delta B_{z}-5A-4a)/8,$ (3)

where $\Delta B_{z}$ is the distance in magnetic field units between the
lowest and the highest transition lines in the ODMR signal.

The above expression enables the determination of the $D$ parameter for
$D\gtrapprox 75$ neV. Below $75$ neV, the transitions between other states are
overlapping, and resolving the lines becomes cumbersome. The difference
between $D$ value calculated from equation 3 for $75$ neV $<D<1000$ neV and
obtained from numerical solving of the spin Hamiltonian for $f=15.6$ GHz is
below $0.05\%$. For higher values of $D$ the difference between the numerical
solution and analytical approximation slowly increases but is still below
$0.1\%$ for $D$ as high as $3000$ neV.

Moreover, the positions of the ODMR signal line crossings in ODMR angular maps
(like those present between B$\approx+20$ mT to B$\approx+30$ mT at
$\theta\approx 30^{\circ}$ and $\theta\approx 50^{\circ}$ in figure 11 or in
figure 11, marked with white arrows) enable for reducing fitting uncertainties
down to a few neV. Similarly, by analyzing analogous features of an in-plane
angular scan we can fit the spin Hamiltonian $E$ parameter which in the case
of QWs should be zero. Indeed within fitting accuracy, the $E$ parameter for
our samples is negligible.

Figure 12 presents the obtained $D$ parameter as a function of deformation. As
it is known Qazzaz  _et al._ (1995) the spin Hamiltonian $D$ parameter is
related to deformation $\varepsilon_{\parallel}$ by formula

$D=-\frac{3}{2}G_{11}\left(1+\frac{2C_{12}}{C_{11}}\right)\varepsilon_{\parallel}$
(4)

where $G_{11}$ is strain spin-lattice coefficient and $C_{11}$, $C_{12}$ are
elastic stiffness constants. The value of $G_{11}$ strain spin-lattice
coefficient obtained from our data is ${G_{11}=(72.2\pm 1.9)}$ neV. This is
significantly larger value than the only previously reported valueCausa and
Tovar (1980) of $G^{\mbox{\footnotesize Causa}}_{11}=(57\pm 1.3)$ neV.

Figure 12: The spin Hamiltonian $D$ parameter obtained from angular ODMR
measurements as a function of deformation for analyzed quantum well samples
(circles). The dashed curve represents the value of G11 for bulk CdMnTe
reported in Ref. Causa and Tovar, 1980. The solid black curve represents the
linear fit to data points measured in this work for Mn2+ in CdTe and
corresponds to ${G_{11}=(72.2\pm 1.9)}$ neV. The dotted lines represent one
standard error deviation.

##  III Summary and conclusions

We have presented a precise method of determination of spin Hamiltonian
parameters by angular resolved ODMR measurements. Due to the high sensitivity
of ODMR measurements (we were probing $\approx$3.8$\times$109 manganese spins
which are almost 3 orders of magnitude less than the limit for standard EPR
setup) we obtained strain-induced axial-symmetry spin Hamiltonian parameter
$D$ for single-layered quantum wells of 10 nm thickness with neV accuracy. By
determination of $D$ values for a series of samples differing by magnesium
content we were able to determine strain spin-lattice coefficient
${G_{11}=(72.2\pm 1.9)}$ neV for Mn2+ ion in CdTe.

##  IV Acknowledgements

This work was supported by the Polish National Science Centre under decisions
DEC-2016/23/B/ST3/03437. The project was carried out with the use of CePT,
CeZaMat, and NLTK infrastructures financed by the European Union - the
European Regional Development Fund within the Operational Programme
”Innovative economy” for 2007 - 2013.

##  Appendix A Temperature dependence of the thermal expansion coefficients of
CdTe and GaAs

The substrate material of analyzed samples is made of GaAs on top of which a
thick CdTe layer is deposited. As the CdTe is deposited at a temperature of
around 560 K, the XRD measurements are performed at 300 K and optical
measurements at 1.6 K it is crucial to take into account the mechanical
properties of the samples. In fact, the presented design is an example of
composite-material behaving similarly to the bimetallic strip, as the thermal
expansion coefficients $\alpha$ of GaAs and CdTe are different. Moreover, the
changes of thermal expansion coefficient of both materials – $\alpha$(T), as
the temperature changes, are different. As a result, accumulated deformation,
due to the cooling, must be evaluated based on full $\alpha$(T) curves. The
values of $\alpha$ are available in the literature for both CdTeNovikova
(1960); Williams  _et al._ (1969); Browder and Ballard (1972); Novikova
(1974); Bagot  _et al._ (1993) and GaAsSparks and Swenson (1967); Feder and
Light (1968, 1968); Novikova (1974); Soma  _et al._ (1982); Sirota  _et al._
(1984). The additional deformation present in CdTe layer can be then
calculated as

$\Delta\varepsilon_{\parallel}^{\textrm{temp.}}=\int_{T_{0}}^{T_{i}}\left(\alpha_{\textrm{GaAs}}(T)-\alpha_{\textrm{CdTe}}(T)\right)dT$
(5)

For convinient integration we interpolated literature data using an empirical
formula:

$\alpha(T)\approx\sum_{i=1}^{i=3}G_{i}(g_{i},s_{i},T)$ (6)

where

$G_{i}(g_{i},s_{i},T)=g_{i}\left(\frac{s_{i}}{T}\right)^{2}\frac{e^{\frac{s_{i}}{T}}}{(e^{\frac{s_{i}}{T}}-1)^{2}}$
(7)

Full theory and formulas are presented elsewhereBagot  _et al._ (1993).
However, only the combination of $\alpha$ curves from ref. Bagot  _et al._ ,
1993 for CdTe with ref. Novikova, 1974 GaAs gives deformation value that
agrees with values obtained from XRD measurements. As the temperature of the
substrate during the MBE growth of the sample was known within 10 K accuracy
we evaluate actual growth temperature $T_{0}$ as a value that results in the
$\varepsilon_{\parallel}=-0.389\text{\textperthousand}$ measured with XRD.
This procedure gives $T_{0}=552$ K.

Figure 13: The upper panel presents linear thermal expansion coefficients of CdTe and GaAs as a function of temperature. Data for CdTe and GaAs were taken from papers of Bagot et al.Bagot  _et al._ (1993) and Novikova et al.Novikova (1974) respectively. Curves are best-fit results of fitting equation 6 to the data. The lower panel shows calculation results of the additional deformation $\Delta\varepsilon_{\parallel}^{\textrm{temp.}}$ present in the QW samples due to the difference in thermal expansion coefficients of CdTe material and GaAs substrate material. The $T_{0}$ is the temperature related to the MBE growth process. Table 2: Numerical values of si, gi coefficients used for reproduction of the linear thermal expansion coefficients $\alpha(T)$ of CdTe and GaAs with eq. 6. material |  | i |  | si (K) | gi (K-1)
---|---|---|---|---|---
CdTe |  | 1 |  | $222.020\,98$ | $883.964\,17$
|  | 2 |  | $-40.900\,23$ | $-384.867\,15$
|  | 3 |  | $-13\,241.61$ | $1403.506\,96$
GaAs |  | 1 |  | $3562.248\,38$ | $640.974\,75$
|  | 2 |  | $-99.546\,97$ | $-44.505\,16$
|  | 3 |  | $371.479\,81$ | $709.906\,72$

##  Appendix B Temperature dependence of the CdTe elastic stiffness constants

R. D. Greenough and S. B. Palmer in ref. Greenough and Palmer, 1973 presented
experimentally determined linear combinations of the CdTe elastic parameters
values from 300 K to 4.2 K. From their data one can calculate $2C_{12}/C_{11}$
ratio. In this paper we evaluate $2C_{12}/C_{11}$ ratio using 2 components of
equation 6 and added constant $C_{0}$. We obtain the best reproduction of data
from Greenough and Palmer, 1973 for g${}_{1}=-0.0092$ K-1, g${}_{2}=-0.015$
K-1, s${}_{1}=-516$ K, s${}_{2}=1415$ K and $C_{0}=1.39886$. For temperatures
below 50 K, as the temperature decreases, the value of $C_{0}$ approaches the
$2C_{12}/C_{11}$ ratio. For 300 K and 1.6 K we find
$2C_{12}/C_{11}|_{300~{}\mathrm{K}}=1.38859$ and
$2C_{12}/C_{11}|_{1.6~{}\mathrm{K}}=1.39886$.

##  Appendix C Numerical simulation of QW levels

Theoretical calculations of energy levels in QWs were performed by solving a
1D Schrodinger equation of carrier confined in the square quantum well. The
effective masses in perpendicular direction are assumed $0.1m_{0}$,
$0.63m_{0}$, and $0.13m_{0}$ for electron, heavy and light hole, respectively
where $m_{0}$ is the free electron mass. The potential was obtained from
composition profiles. The QW layer was composed of
$\mathrm{Cd}_{1-y}\mathrm{Mn}_{y}\mathrm{Te}$, sandwiched between the
$\mathrm{Cd}_{1-x}\mathrm{Mg}{{}_{x}}\mathrm{Te}$ barriers, the chemical
potential $A_{C}$ and $A_{V}$ for conduction and valence band is defined by
the energy gap difference with a relative valence band offset $Q_{V}=0.4$
Wojtowicz  _et al._ (1996).

We used bandgap versus composition relationGaj and Kossut (2010):

$E_{g}^{\mathrm{Cd}_{1-x}\mathrm{Mg}_{x}\mathrm{Te}}(x)=Bx+E_{g}^{\mathrm{CdTe}}=(1.850x+1.606)\mathrm{~{}eV}$
(8)

for the barrier material and for the QW:

$E_{g}^{\mathrm{Cd}_{1-y}\mathrm{Mn}_{y}\mathrm{Te}}(y)=By+E_{g}^{\mathrm{CdTe}}=(1.563y+1.606)\mathrm{~{}eV}.$
(9)

The difference of lattice constants of buffer and barrier layers Dynowska  _et
al._ (1999) leads to the presence of additional deformation potential
$A_{i}^{\mathrm{deform}}=\alpha_{i}\varepsilon_{\parallel}$, where
$i=\mathrm{E,HH,LH}$ denotes electron, heavy and light hole. The $\alpha_{i}$
parameters are corresponding to the $a^{\prime}$ (hydrostatic deformation
potential) and $b$ (shear deformation potential) material constants and
elastic constants ratio. For zinc blend crystals, as CdTe, the modified
potential depths for the electrons, heavy and light holes are Chuang (2009):

$A_{\mathrm{E}}=A_{C}+A_{\mathrm{E}}^{\mathrm{deform}}=A_{C}+\frac{2a^{\prime}}{3}\left(2-\frac{2C_{12}}{C_{11}}\right)\varepsilon_{\parallel},$
(10)

$\displaystyle
A_{\mathrm{HH,LH}}=A_{V}+A_{\mathrm{HH,LH}}^{\mathrm{deform}}=\,\,\,\,$
$\displaystyle=A_{V}+\left[-\frac{a^{\prime}}{3}\left(2-\frac{2C_{12}}{C_{11}}\right)\pm
b\left(1+\frac{2C_{12}}{C_{11}}\right)\right]\varepsilon_{\parallel}.\quad$
(11)

The elastic constants $C_{12}$ and $C_{11}$ are taken from literature
Greenough and Palmer (1973), for this calculation we have used
$2C_{12}/C_{11}=1.39886$, while the parameter $a^{\prime}=-3.85$ eV Peyla  _et
al._ (1992). The light - heavy hole potential difference corresponds only to
the difference in the deformation potential (see equation 1).

The shear deformation potential ${b=(-0.94\pm 0.11)}$ eV was obtained from the
linear fit to the measured light-heavy hole splitting versus calculated
$\varepsilon_{\parallel}$ as it was mentioned before. The exciton binding
energy for heavy hole exciton was set as $18$ meV, while for light hole
exciton $21$ meV, evaluated from the literature Mathieu  _et al._ (1992);
Christol  _et al._ (1993); Gaj  _et al._ (1994). The QW width was equal to
$10.45$ nm, the same for each studied sample. With these assumptions, we have
calculated heavy and light hole exciton energy levels and their energy
difference presented in Fig. 6. Finally, the deformation
$\varepsilon_{\parallel}$ splits heavy and light hole exciton energy levels.

Here we assume that the other parameters as QW width and exciton binding
energy are not vulnerable to the deformation. Moreover, the QW width,
deformation potential corresponding to the $a^{\prime}$ parameter, and the
exciton binding energy affect the exciton energy in a similar way - increasing
or decreasing the state energy. Nonetheless, this additional energy offset
does not affect the heavy-light hole splitting strain dependence which we are
interested in here.

## References

  * Kolm  _et al._ (1957) C. Kolm, S. A. Kulin, and B. L. Averbach, Physical Review 108, 965 (1957).
  * Bassignana  _et al._ (1997) I. C. Bassignana, D. A. Macquistan, G. C. Hillier, R. Streater, D. Beckett, A. Majeed, and C. Miner, Journal of Crystal Growth 178, 445 (1997).
  * Usuda and Ando (1996) K. Usuda and M. Ando, Journal of Applied Physics 80, 1352 (1996).
  * Qazzaz  _et al._ (1995) M. Qazzaz, G. Yang, S. H. Xin, L. Montes, H. Luo, and J. K. Furdyna, Solid State Communications 96, 405 (1995).
  * Causa and Tovar (1980) M. T. Causa and M. Tovar, Physics Letters 77, 473 (1980).
  * Merle d’Aubigné  _et al._ (1990) Y. Merle d’Aubigné, H. Mariette, N. Magnea, H. Tuffigo, R. Cox, G. Lentz, L. Dang, J. Pautrat, and A. Wasiela, Journal of Crystal Growth 101, 650 (1990).
  * Thomas (1961) D. G. Thomas, Journal of Applied Physics 32, 2298 (1961).
  * Mathieu  _et al._ (1988) H. Mathieu, J. Allegre, A. Chatt, P. Lefebvre, and J. P. Faurie, Physical Review B 38, 7740 (1988).
  * Peyla  _et al._ (1992) P. Peyla, Y. Merle d’Aubigné, A. Wasiela, R. Romestain, H. Mariette, M. D. Sturge, N. Magnea, and H. Tuffigo, Physical Review B 46, 1557 (1992).
  * Cherbunin  _et al._ (2020) R. V. Cherbunin, V. M. Litviak, I. I. Ryzhov, A. V. Koudinov, S. Elsässer, A. Knapp, T. Kiessling, J. Geurts, S. Chusnutdinow, T. Wojtowicz, and G. Karczewski, Physical Review B 101, 241301 (2020).
  * Lafuente-Sampietro  _et al._ (2017) A. Lafuente-Sampietro, H. Boukari, and L. Besombes, Physical Review B 95 (2017), 10.1103/PhysRevB.95.245308.
  * Dinu  _et al._ (2017) I. V. Dinu, V. Moldoveanu, R. Dragomir, and B. Tanatar, Physica Status Solidi (B) 254, 1600800 (2017).
  * Moldoveanu  _et al._ (2016) V. Moldoveanu, I. V. Dinu, R. Dragomir, and B. Tanatar, Physical Review B 93, 165421 (2016).
  * Lafuente-Sampietro  _et al._ (2016) A. Lafuente-Sampietro, H. Boukari, and L. Besombes, Physica Status Solidi (C) 13, 546 (2016).
  * Lafuente-Sampietro  _et al._ (2015) A. Lafuente-Sampietro, H. Boukari, and L. Besombes, Physical Review B 92 (2015), 10.1103/PhysRevB.92.081305.
  * Goryca  _et al._ (2015) M. Goryca, M. Koperski, T. Smoleński, Ł. Cywiński, P. Wojnar, P. Plochocka, M. Potemski, and P. Kossacki, Physical Review B 92 (2015), 10.1103/PhysRevB.92.045412.
  * Cronenberger  _et al._ (2015) S. Cronenberger, D. Scalbert, D. Ferrand, H. Boukari, and J. Cibert, Nature Communications 6, 8121 (2015).
  * Kobak  _et al._ (2014) J. Kobak, T. Smoleński, M. Goryca, M. Papaj, K. Gietka, A. Bogucki, M. Koperski, J.-G. Rousset, J. Suffczyński, E. Janik, M. Nawrocki, A. Golnik, P. Kossacki, and W. Pacuski, Nature Communications 5 (2014), 10.1038/ncomms4191.
  * Varghese  _et al._ (2014) B. Varghese, H. Boukari, and L. Besombes, Physical Review B 90 (2014), 10.1103/PhysRevB.90.115307.
  * Besombes and Boukari (2014) L. Besombes and H. Boukari, Physical Review B 89 (2014), 10.1103/PhysRevB.89.085315.
  * Goryca  _et al._ (2014) M. Goryca, M. Koperski, P. Wojnar, T. Smoleński, T. Kazimierczuk, A. Golnik, and P. Kossacki, Physical Review Letters 113, 227202 (2014).
  * Cronenberger  _et al._ (2013) S. Cronenberger, M. Vladimirova, S. V. Andreev, M. B. Lifshits, and D. Scalbert, Physical Review Letters 110, 077403 (2013).
  * Jamet  _et al._ (2013) S. Jamet, H. Boukari, and L. Besombes, Physical Review B 87 (2013), 10.1103/PhysRevB.87.245306.
  * Chen  _et al._ (2013) Y. S. Chen, M. Wiater, G. Karczewski, T. Wojtowicz, and G. Bacher, Physical Review B 87, 155301 (2013).
  * Besombes  _et al._ (2012) L. Besombes, C. L. Cao, S. Jamet, H. Boukari, and J. Fernández-Rossier, Physical Review B 86 (2012), 10.1103/PhysRevB.86.165306.
  * Goryca  _et al._ (2009a) M. Goryca, T. Kazimierczuk, M. Nawrocki, A. Golnik, J. A. Gaj, P. Kossacki, P. Wojnar, and G. Karczewski, Physical Review Letters 103, 087401 (2009a).
  * Goryca  _et al._ (2009b) M. Goryca, D. Ferrand, P. Kossacki, M. Nawrocki, W. Pacuski, W. Maślana, J. A. Gaj, S. Tatarenko, J. Cibert, T. Wojtowicz, and G. Karczewski, Physical Review Letters 102, 046408 (2009b).
  * Le Gall  _et al._ (2009) C. Le Gall, L. Besombes, H. Boukari, R. Kolodka, J. Cibert, and H. Mariette, Physical Review Letters 102 (2009), 10.1103/PhysRevLett.102.127402.
  * Hartmann  _et al._ (1996) J. M. Hartmann, J. Cibert, F. Kany, H. Mariette, M. Charleux, P. Alleysson, R. Langer, and G. Feuillet, Journal of Applied Physics 80, 6257 (1996).
  * LeBlanc  _et al._ (2017) E. G. LeBlanc, M. Edirisooriya, O. S. Ogedengbe, O. C. Noriega, P. A. R. D. Jayathilaka, S. Rab, C. H. Swartz, D. R. Diercks, G. L. Burton, B. P. Gorman, A. Wang, T. M. Barnes, and T. H. Myers, Journal of Electronic Materials 46, 5379 (2017).
  * Cibert  _et al._ (1990) J. Cibert, Y. Gobil, L. S. Dang, S. Tatarenko, G. Feuillet, P. H. Jouneau, and K. Saminadayar, Applied Physics Letters 56, 292 (1990).
  * Cibert  _et al._ (1991) J. Cibert, R. André, C. Deshayes, G. Feuillet, P. H. Jouneau, L. S. Dang, R. Mallard, A. Nahmani, K. Saminadayar, and S. Tatarenko, Superlattices and Microstructures 9, 271 (1991).
  * Gaj  _et al._ (1994) J. A. Gaj, W. Grieshaber, C. Bodin-Deshayes, J. Cibert, G. Feuillet, Y. Merle d’Aubigné, and A. Wasiela, Physical Review B 50, 5512 (1994).
  * Waag  _et al._ (1993) A. Waag, H. Heinke, S. Scholl, C. R. Becker, and G. Landwehr, Journal of Crystal Growth 131, 607 (1993).
  * Kossacki  _et al._ (1997) P. Kossacki, N. The Khoi, J. A. Gaj, M. Kutrowski, M. Surma, G. Karczewski, T. Wojtowicz, and M. Guziewicz, Solid State Communications 103, 619 (1997).
  * Kossacki (2003) P. Kossacki, Journal of Physics: Condensed Matter 15, R471 (2003).
  * Kossacki  _et al._ (2004) P. Kossacki, H. Boukari, M. Bertolini, D. Ferrand, J. Cibert, S. Tatarenko, J. A. Gaj, B. Deveaud, V. Ciulin, and M. Potemski, Physical Review B 70 (2004), 10.1103/PhysRevB.70.195337.
  * Allègre  _et al._ (1990) J. Allègre, B. Gil, J. Calatayud, and H. Mathieu, Journal of Crystal Growth 101, 603 (1990).
  * Boley  _et al._ (1993) M. S. Boley, R. J. Thomas, M. Chandrasekhar, H. R. Chandrasekhar, A. K. Ramdas, M. Kobayashi, and R. L. Gunshor, Journal of Applied Physics 74, 4136 (1993).
  * Bhattacharjee and Fishman (1998) A. K. Bhattacharjee and G. Fishman, Journal of Crystal Growth 184-185, 942 (1998).
  * Horning and Staudenmann (1987) R. D. Horning and J.-L. Staudenmann, Applied Physics Letters 50, 1482 (1987).
  * Tatsuoka  _et al._ (1989) H. Tatsuoka, H. Kuwabara, H. Fujiyasu, and Y. Nakanishi, Journal of Applied Physics 65, 2073 (1989).
  * Tatsuoka  _et al._ (1990) H. Tatsuoka, H. Kuwabara, Y. Nakanishi, and H. Fujiyasu, Journal of Applied Physics 67, 6860 (1990).
  * Tatsuoka  _et al._ (1991) H. Tatsuoka, H. Kuwabara, Y. Nakanishi, and H. Fujiyasu, Thin Solid Films 201, 59 (1991).
  * Chuang (2009) S. L. Chuang,  _Physics of Photonic Devices_ , 2nd ed., Wiley Series in Pure and Applied Optics (John Wiley & Sons, Hoboken, N.J, 2009).
  * van Gisbergen  _et al._ (1993) S. J. C. H. M. van Gisbergen, M. Godlewski, R. R. Gal/a̧zka, T. Gregorkiewicz, C. A. J. Ammerlaan, and N. T. Khoi, Physical Review B 48, 11767 (1993).
  * Ivanov  _et al._ (2001) V. Ivanov, M. Godlewski, S. Yatsunenko, A. Khachapuridze, M. Li, and Z. Gołacki, Acta Physica Polonica A 100, 351 (2001).
  * Godlewski  _et al._ (2002) M. Godlewski, V. Y. Ivanov, A. Khachapuridze, and S. Yatsunenko, physica status solidi (b) 229, 533 (2002).
  * Khachapuridze  _et al._ (2002) A. Khachapuridze, V. Ivanov, M. Godlewski, and G. Karczewski, Acta Physica Polonica A 102, 643 (2002).
  * Godlewski  _et al._ (2004) M. Godlewski, S. Yatsunenko, A. Khachapuridze, V. Y. Ivanov, Z. Gołacki, G. Karczewski, P. J. Bergman, P. J. Klar, W. Heimbrodt, and M. R. Phillips, Journal of Alloys and Compounds Proceedings of the 4th International Spring Workshop on Spectroscopy, Structure and Synthesis of Rare Earth Systems, 380, 45 (2004).
  * Smith  _et al._ (2005) L. C. Smith, S. J. Bingham, J. J. Davies, and D. Wolverson, Applied Physics Letters 87, 202101 (2005).
  * Godlewski  _et al._ (2008) M. Godlewski, V. Y. Ivanov, and S. Yatsunenko, Optical Materials Proceedings of the International Workshop on Advanced Spectroscopy and Optical Materials, 30, 669 (2008).
  * Baranov  _et al._ (2009) P. G. Baranov, N. G. Romanov, D. O. Tolmachev, R. A. Babunts, B. R. Namozov, Y. G. Kusrayev, I. V. Sedova, S. V. Sorokin, and S. V. Ivanov, JETP Letters 88, 631 (2009).
  * Debus  _et al._ (2016) J. Debus, V. Y. Ivanov, S. M. Ryabchenko, D. R. Yakovlev, A. A. Maksimov, Y. G. Semenov, D. Braukmann, J. Rautert, U. Löw, M. Godlewski, A. Waag, and M. Bayer, Physical Review B 93, 195307 (2016).
  * O. Tolmachev  _et al._ (2020) D. O. Tolmachev, V. Yu. Ivanov, D. R. Yakovlev, E. V. Shornikova, B. Witkowski, S. Shendre, F. Isik, S. Delikani, H. Volkan Demir, and M. Bayer, Nanoscale 12, 21932 (2020).
  * Shornikova  _et al._ (2020) E. V. Shornikova, D. R. Yakovlev, D. O. Tolmachev, V. Y. Ivanov, I. V. Kalitukha, V. F. Sapega, D. Kudlacik, Y. G. Kusrayev, A. A. Golovatenko, S. Shendre, S. Delikanli, H. V. Demir, and M. Bayer, ACS Nano 14, 9032 (2020).
  * Komarov  _et al._ (1977) A. V. Komarov, S. Ryabchenko, O. Terletskii, and I. Zheru, Journal of Experimantal and Theoretical Physics 46, 318 (1977).
  * Ivanov  _et al._ (2008) V. Y. Ivanov, M. Godlewski, D. R. Yakovlev, M. K. Kneip, M. Bayer, S. M. Ryabchenko, and A. Waag, Physical Review B 78 (2008), 10.1103/PhysRevB.78.085322.
  * Ivanov and Godlewski (2010) V. Y. Ivanov and M. Godlewski, Applied Magnetic Resonance 39, 31 (2010).
  * Sadowski  _et al._ (2003) M. L. Sadowski, M. Byszewski, M. Potemski, A. Sachrajda, and G. Karczewski, Applied Physics Letters 82, 3719 (2003).
  * Byszewski  _et al._ (2004) M. Byszewski, D. Plantier, M. Sadowski, M. Potemski, A. Sachrajda, Z. Wilamowski, and G. Karczewski, Physica E: Low-dimensional Systems and Nanostructures 22, 652 (2004).
  * Tolmachev  _et al._ (2010) D. O. Tolmachev, R. A. Babunts, N. G. Romanov, P. G. Baranov, B. R. Namozov, Y. G. Kusrayev, S. Lee, M. Dobrowolska, and J. K. Furdyna, Physica Status Solidi (B) 247, 1511 (2010).
  * Tolmachev  _et al._ (2012) D. O. Tolmachev, A. S. Gurin, N. G. Romanov, A. G. Badalyan, R. A. Babunts, P. G. Baranov, B. R. Namozov, and Y. G. Kusrayev, JETP Letters 96, 231 (2012).
  * Gurin  _et al._ (2015) A. S. Gurin, D. O. Tolmachev, N. G. Romanov, B. R. Namozov, P. G. Baranov, Y. G. Kusrayev, and G. Karczewski, JETP Letters 102, 230 (2015).
  * Keller  _et al._ (2001) D. Keller, D. R. Yakovlev, B. König, W. Ossau, T. Gruber, A. Waag, L. W. Molenkamp, and A. V. Scherbakov, Physical Review B 65, 035313 (2001).
  * Seufert  _et al._ (2001) J. Seufert, G. Bacher, M. Scheibner, A. Forchel, S. Lee, M. Dobrowolska, and J. K. Furdyna, Physical Review Letters 88, 027402 (2001).
  * Yakovlev  _et al._ (2004) D. R. Yakovlev, M. Kneip, A. A. Maksimov, I. I. Tartakovskii, M. Bayer, D. Keller, W. Ossau, L. W. Molenkamp, A. V. Scherbakov, A. V. Akimov, and A. Waag, Physica Status Solidi (C) 1, 989 (2004).
  * Hundt  _et al._ (2005) A. Hundt, J. Puls, A. V. Akimov, Y. H. Fan, and F. Henneberger, Physical Review B 72, 033304 (2005).
  * Abragam and Bleaney (2012) A. Abragam and B. Bleaney,  _Electron Paramagnetic Resonance of Transition Ions_ , Oxford Classic Texts in the Physical Sciences (Oxford University Press, Oxford, 2012).
  * Novikova (1960) S. Novikova, Fizika Tverdogo Tela 2, 2341 (1960).
  * Williams  _et al._ (1969) M. Williams, R. Tomlinson, and M. Hampshire, Solid State Communications 7, 1831 (1969).
  * Browder and Ballard (1972) J. S. Browder and S. S. Ballard, Applied Optics 11, 841 (1972).
  * Novikova (1974) S. Novikova,  _Teplovoe Rasshirenie Tverdyh Tel_ (Academy of Sciences of the Soviet Union, IzdatelĘšstvo "Nauka", Moscow, 1974).
  * Bagot  _et al._ (1993) D. Bagot, R. Granger, and S. Rolland, Physica Status Solidi (B) 177, 295 (1993).
  * Sparks and Swenson (1967) P. W. Sparks and C. A. Swenson, Physical Review 163, 779 (1967).
  * Feder and Light (1968) R. Feder and T. Light, Journal of Applied Physics 39, 4870 (1968).
  * Soma  _et al._ (1982) T. Soma, J. Satoh, and H. Matsuo, Solid State Communications 42, 889 (1982).
  * Sirota  _et al._ (1984) N. N. Sirota, A. M. Antyukhov, and A. A. Sidorov, Doklady Akademii Nauk SSSR 277, 1379 (1984).
  * Greenough and Palmer (1973) R. D. Greenough and S. B. Palmer, Journal of Physics D: Applied Physics 6, 587 (1973).
  * Wojtowicz  _et al._ (1996) T. Wojtowicz, M. Kutrowski, M. Surma, K. Kopalko, G. Karczewski, J. Kossut, M. Godlewski, P. Kossacki, and N. The Khoi, Applied Physics Letters 68, 3326 (1996).
  * Gaj and Kossut (2010) J. A. Gaj and J. Kossut, eds.,  _Introduction to the Physics of Diluted Magnetic Semiconductors_, Springer Series in Materials Science (Springer-Verlag, Berlin Heidelberg, 2010).
  * Dynowska  _et al._ (1999) E. Dynowska, E. Janik, J. Ba̧k-Misiuk, J. Domagała, T. Wojtowicz, and J. Kossut, Journal of Alloys and Compounds 286, 276 (1999).
  * Mathieu  _et al._ (1992) H. Mathieu, P. Lefebvre, and P. Christol, Physical Review B 46, 4092 (1992).
  * Christol  _et al._ (1993) P. Christol, P. Lefebvre, and H. Mathieu, Journal of Applied Physics 74, 5626 (1993).


the Brownian increment. Hence, the discretization factor $\Delta t$ only
appears in the deterministic part.

The following result follows directly from Theorem 6.12 and [22, Theorem
1.46].

######  Corollary 6.14.

Suppose that all assumptions of Theorem 6.12 hold. We assume that
$\bar{u}\in{{\mathscr{U}}_{\mathsf{ad}}}$ is a local solution to (29). Then,
first-order necessary optimality conditions are given by the variational
inequality

$\displaystyle\bigg{\langle}\kappa\bar{u}+{\mathbb{E}}\Big{[}-a_{u}({\bm{\bar{X}}},\bar{u},\cdot)^{\top}\bm{\bar{\Lambda}}-\sum_{j=1}^{m}b_{ju}({\bm{\bar{X}}},\bar{u},\cdot)^{\top}\bm{\bar{Z}}_{j}\Big{]},u-\bar{u}\bigg{\rangle}_{{\mathscr{U}}}\geq
0\quad\text{for all }u\in{{\mathscr{U}}_{\mathsf{ad}}},$

where ${\bm{\bar{X}}}={\mathcal{S}}(\bar{u})$ holds and
$({\bm{\Lambda}},{\bm{Z}})$ solves (91) for
$({\bm{X}},u)=({\bm{\bar{X}}},\bar{u})$.

##  7 Numerical experiments

In this section, we explain our optimization strategy. Afterwards, we perform
numerical experiments. In the first one, we find an optimal control such that
the mean and variance of the governing process follow precisely the prescribed
data. The governing process in this case is the well-known mean-reverting
Ornstein-Uhlenbeck process. In the second one, we calibrate a model consisting
of systems of SDEs by finding optimal parameters. The specific model is the
so-called Stochastic Prandl-Tomlinson (SPT) model that is used to study
microrheological processes of viscous fluids.

Let us mention that in both examples our cost functionals are not of the
specific form introduced in Example 4.2.

#### Optimization strategy

Exploiting the discrete optimality system and using the Euler-Maruyama (EM)
scheme [21, 36], we construct a stochastic gradient method that is summarized
in Algorithm 1. During the process of solving the model and adjoint SDE, we
utilize the parallelization technique of multi-threading in order to speed up
the calculations. This is possible since each realization is independent of
the other ones. While updating the control, we apply a projection onto
${{\mathscr{U}}_{\mathsf{ad}}}$ or
${{\mathscr{U}}^{\mathsf{c}}_{\mathsf{ad}}}$ utilizing the standard projection
operator $\mathcal{P}_{{{\mathscr{U}}_{\mathsf{ad}}}}$ that is based on a
componentwise projection.

Hence, we apply a stochastic gradient method using fixed batch sizes
$M^{\ell}$. Convergence properties can be found in [6]. In our examples, we
use with some initial $s_{0}>0$ the stepsize rule

(96) $\displaystyle s_{\ell}=s_{0}/\ell\quad\text{with }\ell\in{\mathbb{N}}$

for the stochastic gradient method.

Algorithm 1 Stochastic Gradient Method for Stochastic Optimal Control

0: Initial discrete control guess $u^{0}=(u^{0}_{1},\ldots,u^{0}_{r})^{\top}$
in case of ${\mathscr{U}}=[{u_{\mathsf{a}}},{u_{\mathsf{b}}}]$ (calibration
problem) or ${\bm{U}}^{h,0}=[\mathrm{u}_{0}^{0}|\ldots|\mathrm{u}_{N-1}^{0}]$
in case of ${\mathscr{U}}=L^{2}(0,T;{\mathbb{R}}^{r})$ (time-dependent
controls), tolerance $tol>0$, maximum iteration depth
$\ell_{\max}\in{\mathbb{N}}$;

1: Set $\ell=0$, initialize $E\gg tol$;

2: while $E>tol$ and $\ell<\ell_{\max}$ do

3: Generate $M^{\ell}\in{\mathbb{N}}$ random numbers (call the set of the
random numbers ${\mathscr{M}}^{\ell}\in{\mathbb{N}}^{M^{\ell}}$);

4: Solve the discretized state model using the EM scheme;

5: Solve the associated adjoint model using the EM scheme;

6: Assemble the gradient $\nabla\hat{J}^{h}(u^{\ell})$ and
$\nabla\hat{J}^{h}({\bm{U}}^{h,\ell})$, respectively, of the reduced cost and
compute $E=\|\nabla\hat{J}^{h}(u^{\ell+1})\|_{{\mathscr{U}}^{\mathsf{c}}}$ and
$E=\|\nabla\hat{J}^{h}({\bm{U}}^{h,\ell+1})\|_{{\mathscr{U}}^{h}}$,
respectively;

7: Calculate stepsize $s_{\ell}$ using (96);

8: Determine a new control $u^{\ell+1}$ and ${\bm{U}}^{h,\ell+1}$,
respectively, by applying a projected stochastic gradient step;

9: Set $\ell=\ell+1$;

10: end while

11: return $u^{\ell}$ and ${\bm{U}}^{h,\ell}$, respectively.

######  Remark 7.1.

The random numbers ${\mathscr{M}}^{\ell}$ in Algorithm 1 change in every
optimization iteration, but stay the same within an iteration.

###  7.1 Time-dependent Ornstein-Uhlenbeck process

In the context of (14) we have $d=1$, $m=1$, $r=2$ and therefore
${\mathscr{U}}\colonequals L^{2}(0,T;{\mathbb{R}}^{2})$. For a fixed parameter
$\theta>0$, the stochastic process $X(t)\in L^{2}(\Omega;{\mathbb{R}})$ solves
the linear SDE

(97)
$\displaystyle\mathrm{d}X(t)=\theta\big{(}u_{1}(t)-X(t)\big{)}\,\mathrm{d}t+u_{2}(t)\,\mathrm{d}B(t)\text{
for }t\in(0,T],\quad X(0)=X_{\circ}$

with an initial condition $X_{\circ}\in L^{2}(\Omega)$ and the time-dependent
coefficients $u=(u_{1},u_{2})\in{\mathscr{U}}$. We set

$\displaystyle a(x,u,t)=\theta\big{(}u_{1}-x\big{)},\quad
b(x,u,t)=u_{2}\quad\text{for
}(x,u,t)\in{\mathbb{R}}\times{\mathbb{R}}^{r}\times[0,T]\text{ and
}u=(u_{1},u_{2}).$

Then, (97) can be expressed in the form (14). Recall that for $X\in
L^{2}(\Omega)$ the variance is defined as

$\displaystyle{\mathbb{V}}[X]\colonequals{\mathbb{E}}\big{[}X-{\mathbb{E}}[X]^{2}\big{]}=\int_{\Omega}\big{(}X-{\mathbb{E}}[X]\big{)}^{2}\,\mathrm{d}{\mathbb{P}}=\int_{\Omega}\big{(}X(\omega)-{\mathbb{E}}[X]\big{)}^{2}\,\mathrm{d}{\mathbb{P}}(\omega).$

Now we introduce the cost functional as

(98) $\displaystyle\begin{aligned}
J(X,u)&\colonequals\frac{1}{2}\,\big{\|}{\mathbb{E}}[X(\cdot)]-{\eta^{\mathsf{d}}}\big{\|}_{L^{2}(0,T)}^{2}+\frac{1}{2}\,\big{\|}\mathbb{V}[X(\cdot)]-{\sigma^{\mathsf{d}}}\big{\|}_{L^{2}(0,T)}^{2}\\\
&\quad~{}+\frac{1}{2}\,\big{|}{\mathbb{E}}[X(T)]-{\eta^{\mathsf{d}}_{T}}\big{|}^{2}+\frac{1}{2}\,\big{|}\mathbb{V}[X(T)]-{\sigma^{\mathsf{d}}_{T}}\big{|}^{2}+\frac{\kappa}{2}\,{\|u\|}_{\mathscr{U}}^{2}\end{aligned}$

with

$\displaystyle{\eta^{\mathsf{d}}}(t)=\sin\Big{(}\frac{2\pi
t}{T}\Big{)}-1,\quad{\sigma^{\mathsf{d}}}(t)=0.2\bigg{(}\cos\Big{(}\frac{2\pi
t}{T}\Big{)}+2\bigg{)},\quad{\eta^{\mathsf{d}}_{T}}={\eta^{\mathsf{d}}}(T),\quad{\sigma^{\mathsf{d}}_{T}}={\sigma^{\mathsf{d}}}(T)$

for $t\in[0,T]$. Setting for $t\in[0,T]$ and $X\in{\mathscr{X}}$

$\displaystyle j(X)$
$\displaystyle=\frac{1}{2}\,\big{\|}{\mathbb{E}}[X(\cdot)]-{\eta^{\mathsf{d}}}\big{\|}_{L^{2}(0,T)}^{2}+\frac{1}{2}\,\big{\|}{\mathbb{V}}[X(\cdot)]-{\sigma^{\mathsf{d}}}\big{\|}^{2}_{L^{2}(0,T)},$
$\displaystyle j_{T}(X)$
$\displaystyle=\frac{1}{2}\,\big{\|}{\mathbb{E}}[X(T)]-{\eta^{\mathsf{d}}_{T}}\big{\|}_{L^{2}(0,T)}^{2}+\frac{1}{2}\,\big{\|}{\mathbb{V}}[X(T)]-{\sigma^{\mathsf{d}}_{T}}\big{\|}^{2}_{L^{2}(0,T)}$

the cost (98) and our optimization problem can be expressed as (25) and (27),
respectively.

######  Remark 7.2.

  1. 1)

Note that the coefficient functions $a$ and $b$ are independent of $t$ and
affine linear in $x$ and $u$. In particular, both coefficients are Lipschitz
continuous in $x$ and $u$. Moreover, $a$ and $b$ are time-independent.
Furthermore, (13) holds for $u\in{\mathscr{U}}$. Hence, Assumptions 3.2 and
4.3 are fulfilled.

  2. 2)

For the cost functional it holds that it is continuous as a composition of
continuous functions and moreover bounded from below by zero. Furthermore, the
parts of the functional only containing the expected value are convex as a
composition of non-decreasing functions and a convex function. However, since
the variance is not a convex function, the parts with the variance are not
convex. Consequently, Assumption 4.1 is only partially satisfied.

  3. 3)

Note that $a$ and $b$ are continuously differentiable and hence Assumption
5.1-2) is fulfilled. Furthermore, $a_{x}(x,u(t),t)=-\theta$ and
$b_{x}(x,u(t),t)=0$ hold for $(x,u)\in{\mathbb{R}}\times{\mathscr{U}}$ and
$t\in[0,T]$. Thus, (70) is clearly valid by choosing
$L_{\mathfrak{ab}}=\theta$ and $L_{\mathfrak{ab}}=0$, respectively. Thus,
Assumption 6.4 is fulfilled.

Utilizing the notation introduced in Section 5.3, it turns out that for
$\mu=1,\ldots,M$ and $\nu=0,\ldots,N-1$, the first-order necessary optimality
system of the discrete optimization problem is given by

(99a) $\displaystyle\left\\{\begin{aligned}
X^{\mu}_{\nu+1}&=X^{\mu}_{\nu}+\big{(}\theta({\bm{U}}^{h}_{1\nu}-X^{\mu}_{\nu})\big{)}\,\Delta
t+{\bm{U}}^{h}_{2\nu}\,\Delta B^{\mu}_{\nu},\\\
X^{\mu}_{0}&=X_{\circ}^{\mu},\end{aligned}\right.$ (99b)
$\displaystyle\left\\{\begin{aligned}
\Lambda^{\mu}_{\nu}&=\Lambda_{\nu+1}^{\mu}-\theta\Lambda_{\nu+1}^{\mu}\,\Delta
t+\Delta
t\left({\mathbb{E}}^{M}\big{[}X_{\nu}^{M}\big{]}-\eta^{\mathrm{d}}_{\nu}\right)\\\
&\quad+2\Delta
t\left({\mathbb{V}}^{M}\big{[}X_{\nu}^{M}\big{]}-\sigma_{\nu}^{\mathsf{d}}\right)\left(X^{\mu}_{\nu}-{\mathbb{E}}^{M}\big{[}X_{\nu}^{M}\big{]}\right),\\\
\Lambda^{\mu}_{N}&=\left({\mathbb{E}}^{M}\big{[}X_{N}^{h}\big{]}-{\eta^{\mathsf{d}}_{T}}\right)+\Delta
t\left({\mathbb{E}}^{M}\big{[}X_{N}^{h}\big{]}-\mu^{\mathrm{d}}_{N}\right)\\\
&\quad+2\Delta
t\left({\mathbb{V}}^{M}\big{[}X_{N}^{h}\big{]}-\sigma_{N}^{\mathsf{d}}\right)\left(X^{\mu}_{N}-{\mathbb{E}}^{M}\big{[}X_{N}^{h}\big{]}\right)\\\
&\quad+2\left({\mathbb{V}}^{M}\big{[}X_{N}^{h}\big{]}-\sigma_{N}^{\mathsf{d}}\right)\left(X^{\mu}_{N}-{\mathbb{E}}^{M}\big{[}X_{N}^{h}\big{]}\right),\end{aligned}\right.$
(99c)
$\displaystyle\nabla\hat{J}^{h}_{\nu}(\mathrm{u})\colonequals\kappa\,\mathrm{u}_{\nu}-\frac{1}{M}\sum_{\mu=1}^{M}\left(\binom{\theta}{0}\Delta
t+\binom{0}{1}\Delta B^{\mu}_{\nu}\right)\Lambda_{\nu+1}^{\mu}\quad\text{for
}\nu=0,\ldots,N-1$

with

$\displaystyle X^{h}$
$\displaystyle=\big{\\{}X_{\nu}^{\mu}\in{\mathbb{R}}\,\big{|}\,\mu=1,\ldots,M,~{}\nu=0,\ldots,N\big{\\}},$
$\displaystyle{\bm{U}}^{h}$
$\displaystyle=\big{[}\mathrm{u}_{0}|\ldots|\mathrm{u}_{N-1}\big{]}\in{\mathbb{R}}^{2\times
N},\quad\mathrm{u}_{\nu}=\big{(}{\bm{U}}^{h}_{1\nu},{\bm{U}}^{h}_{2\nu}\big{)}^{\top}\in{\mathbb{R}}^{2}~{}(\nu=0,\ldots,N-1),$
$\displaystyle\Lambda^{h}$
$\displaystyle=\big{\\{}\Lambda_{\nu}^{\mu}\in{\mathbb{R}}\,\big{|}\,\mu=1,\ldots,M,~{}\nu=0,\ldots,N\big{\\}}.$

In the sequel, we set the regularization parameter to zero, i.e. $\kappa=0$.
Notice that (99c) can in this case be written as

(100)
$\displaystyle\nabla\hat{J}^{h}_{\nu}(\mathrm{u})={\mathbb{E}}^{M}\left[\binom{\theta\Delta
t\,\Lambda_{\nu+1}^{M}}{\Delta
B_{\nu}\,\Lambda^{M}_{\nu+1}}\right]=\binom{\theta\Delta
t\,{\mathbb{E}}^{M}[\Lambda_{\nu+1}^{M}]}{{\mathbb{E}}^{M}[\Delta
B_{\nu}\Lambda^{M}_{\nu+1}]}.$

Recall that ${\mathbb{E}}[\Delta B_{\nu+1}]=0$. However, since we used the
$\Delta B^{\mu}_{\nu}$ to derive $X_{\nu+1}$ and with this $\Lambda_{\nu+1}$,
it holds that $\Lambda_{\nu+1}$ is not independent of $\Delta B_{\nu}$.
Therefore, the second component of the gradient in (100) is non-zero.

#### Numerical results

In the next figures, we present the results of our optimization strategy for
the Ornstein-Uhlenbeck process (97). For this example, we choose the initial
condition $X_{\circ}$ to obey the normal distribution with mean
$\eta^{\mathsf{d}}(0)$ and standard deviation $\sigma^{\mathsf{d}}(0)$.

Since (97) is a linear SDE, the trajectory of mean and variance can be
calculated as (see, e.g., [35, Chapter 3, Theorem 3.2 and Example 5.2])

(101) $\displaystyle\begin{split}{\mathbb{E}}\big{[}X(t)\big{]}&=e^{-\theta
t}{\mathbb{E}}[X_{\circ}]+e^{-\theta t}\theta\int_{0}^{s}u_{1}(s)e^{\theta
s}\,d\mathrm{d}s,\\\ \mathbb{V}\big{[}X(t)\big{]}&=e^{-2\theta
t}\mathbb{V}\big{[}X_{\circ}\big{]}+e^{-2\theta
t}\int_{0}^{t}u_{2}^{2}(s)e^{2\theta s}\,\mathrm{d}s.\end{split}$

Hence, we can derive the best controls $u$ in the sense that they will lead to
a perfect tracking of the desired trajectories of mean and variance:

(102) $\displaystyle
u_{1}^{*}(t)=(\eta^{\mathsf{d}})^{\prime}(t)\frac{1}{\theta}+\eta^{\mathsf{d}}(t)\quad\text{and}\quad(u_{2}^{*})^{2}(t)=2\theta(\sigma^{\mathsf{d}})(t)+\sqrt{\sigma^{\mathsf{d}}(t)}(\sqrt{\sigma^{\mathsf{d}}})^{\prime}(t).$

The success of the method given in Algorithm 1 is evident in Fig. 1. The
method manages to find controls $u_{1}$ and $u_{2}$ such that the trajectory
of the mean value and the variance of the solution of the model equation
follows the desired one. In Fig. 1(a), the trajectories of all $M$ trials are
plotted together with the calculated and desired moments. In Fig. 1(b), the
control that is calculated with our optimization method is plotted together
with the theoretical control $(u_{1}^{*}(t),u_{2}^{*}(t))$ that leads to a
perfect tracking. One obtains a very good agreement for $u_{1}(t)$. For
$u_{2}(t)$ there is higher noise. However, for more iterations of the
stochastic gradient method, the noise gets smaller. Finally, in Fig. 1(c), we
plot the convergence history of the (relative) functional and the relative
norm of the gradient over the optimization iterations $\ell$. and the
analogous holds for the relative gradient.

(a)

(b)

(c)

Figure 1:  Results for the example of Section 7.1. (a) Trajectories of all
trials (grayed out), desired moments (red) calculated moments (black) while
applying the control calculated by Algorithm 1; (b) Calculated controls $u$
and $u$ together with controls $u_{1}^{*}$ and $u_{2}^{*}$ defined in (102)
for a perfect tracking; (c) Convergence history of relative functional (solid)
and Euclidean norm of the relative gradient (dashed) over the optimization
iterations (cf. (103)).

More specifically, the relative functional is given by

(103)
$\displaystyle\hat{J}^{h}({\bm{U}}^{h,\ell})/|\hat{J}^{h}({\bm{U}}^{h,0})|,$

###  7.2 Stochastic Prantl-Tomlinson (SPT) model

This model is a non-equilibrium bath model in which a tracer particle, also
called colloidal, is assumed to be immersed in a suspension of micelles; see,
e.g., [24, 25, 38] and the references therein. In this setting, the internal
forces between the particles lead to non-Markovian behavior of the movement of
the colloidal due to the memory of the system. The goal is to study the
properties of the suspension by tracing the colloidal.

To model the complex fluid, the colloidal particle is investigated with a
coupling to one or more bath particles that describe the background. The
system of differential equations describing the position of the tracer and $K$
bath particles is given by

(104) $\displaystyle\begin{aligned}
\mathrm{d}X_{1}(t)&=\frac{1}{\gamma_{1}}\bigg{(}-\partial_{x}V_{\mathrm{ext}}(X_{1}(t),t)-\sum_{k=2}^{K+1}V_{\mathrm{int},k}^{\prime}\big{(}X_{1}(t)-X_{k}(t)\big{)}\bigg{)}\,\mathrm{d}t+\frac{1}{\gamma_{1}}\,\mathrm{d}B_{1}(t),\\\
\mathrm{d}X_{k}(t)&=\frac{1}{\gamma_{k}}\,V_{\mathrm{int},k}^{\prime}\big{(}X_{1}(t)-X_{k}(t)\big{)}+\frac{1}{\gamma_{k}}\,\mathrm{d}B_{k}(t)\quad\text{for
}k=2,\ldots,K+1,\end{aligned}$

for $t\in(0,T]$, where for $(x,t)\in\mathbb{R}\times[0,T]$

$\displaystyle
V_{\mathrm{ext}}(x,t)=\frac{\kappa_{\mathrm{ext}}}{2}\,\big{(}x-v_{0}t\big{)}^{2}\qquad\text{and}\qquad
V_{\mathrm{int},k}(x)=V_{0,k}\cos\bigg{(}\frac{2\pi x}{d_{k}}\bigg{)}\text{
for }k=2,\ldots,K+1.$

The system (104) is completed with the initial condition
$X_{i}(0)=X^{eq}_{i,\circ}$, $i=1,\ldots,K+1$, where $X^{eq}_{i,\circ}\in
L^{2}(\Omega;{\mathbb{R}})$ is the equilibrium distribution of the
corresponding particle given the external potential $V_{\mathrm{ext}}(X,0)$.
In system (104), the variable $X_{1}$ denotes the position of the tracer
particle and the $X_{k}$’s the position of the bath particle $p_{k}$,
$k=2,\ldots,K+1$. The $\gamma_{i}>0$, $i=1,\ldots,K+1$, are the friction
coefficients of the corresponding particle and are assumed to be known. The
external potential $V_{\mathrm{ext}}$ may be applied using optical traps
(laser beams) of strength $\kappa_{\mathrm{ext}}>0$ to realize harmonic
potentials. This potential has a trapping effect, which means that the
particle cannot leave the potential easily. This trap is then pulled with a
constant velocity $v_{0}$ through the surrounding medium. The internal forces
between the tracer and the bath particles are denoted by $V_{\mathrm{int},k}$
for $k=2,\ldots,K+1$. The force

$\displaystyle\text{``$\xi_{k}(t)=\dot{B}_{k}(t)$''}\quad\text{for
}t\in[0,T]\text{ and }k=1,\ldots,K+1.$

resembles a random force for which the following mean value relation and
fluctuation-dissipation theorem must hold

(105)
$\displaystyle{\mathbb{E}}\big{[}\xi_{i}(t)\big{]}=0,\quad{\mathbb{E}}\big{[}\xi_{i}(t)\xi_{j}(t^{\prime})\big{]}=2k_{B}\mathrm{T}\gamma_{i}\delta_{ij}\delta(t-t^{\prime})\quad\text{for
}i,j=1,\ldots,K+1,$

where $k_{B}>0$ is Boltzmann’s constant and $\mathrm{T}>0$ the temperature of
the bath. The requirements (105) are fulfilled by Brownian motion, hence we
can apply our setting above.

######  Remark 7.3.

Notice that in [24, 25], the authors did not specify any initial condition
within the SPT model (104). The reason for this is that the system needs to be
in equilibrium. The realize this numerically, we need to integrate over the
transient phase in order to generate sensible data that will be used to
calculate the simulation data. To generate ${\bm{X}}^{eq}_{\circ}$, any
initial condition can be used and the model needs to evolve for a certain
equilibration time $t_{eq}>0$. Then, we use the state of the system after
$t_{eq}$ as ${\bm{X}}^{eq}_{\circ}$. In order to account for a possible high
equilibration time, one needs to resort to parallelization techniques such as
multi-threading that we use in our simulations.

The goal is to identify the pairs
$\\{(V_{0,k},d_{k})\\}_{k=2}^{K+1}\subset{\mathbb{R}}^{2}$ for each bath
particle $p_{k}$, $k=2,\ldots,K+1$, within (104), such that experimental data
is matched as good as possible. For $r=2K$ we define the parameter vector
$u=[u_{1},\ldots,u_{2K}]^{\top}\in{{\mathscr{U}}^{\mathsf{c}}}\colonequals{\mathbb{R}}^{r}$
with

$\displaystyle u_{2k-3}=V_{0,k}\quad\text{and}\quad
u_{2k-2}=\frac{1}{d_{k}}\quad\text{for }k=2,\ldots,K+1.$

It follows that

$\displaystyle V_{\mathrm{int},k}(x)$ $\displaystyle=u_{2k-3}\cos\big{(}2\pi
u_{2k-2}x\big{)}$ $\displaystyle\text{for }k=2,\ldots,K+1,$ $\displaystyle
V_{\mathrm{int},k}^{\prime}(x)$ $\displaystyle=-2\pi
u_{2k-3}u_{2k-2}\sin\big{(}2\pi u_{2k-2}x\big{)}$ $\displaystyle\text{for
}k=2,\ldots,K+1.$

We suppose that ${u_{\mathsf{a}}}\in{{\mathscr{U}}^{\mathsf{c}}}$ satisfies
${u_{\mathsf{a}}}>0$ component-wise in ${\mathbb{R}}^{2K}$. Further, the state
variable consists of the tracer particle $X_{1}$ and the bath particles
$X_{2},\ldots,X_{K+1}$. Hence, for $d=K+1$ and $t\in[0,T]$ we define the
following state vector

$\displaystyle{\bm{X}}(t)\colonequals\big{(}X_{1}(t),\ldots,X_{d}(t)\big{)}^{\top}\in{\mathbb{R}}^{d}.$

The SDE for ${\bm{X}}$ is given by (14), where the coefficient functions $a$
and $b$ are given for every
$(x,u,t)\in{\mathbb{R}}^{d}\times{\mathbb{R}}^{r}\times[0,T]$ as

(106a) $\displaystyle\begin{aligned}
a_{1}(x,u,t)&=\frac{1}{\gamma_{1}}\bigg{(}-\partial_{x}V_{\mathrm{ext}}(x_{1},t)-\sum_{k=2}^{d}V_{\mathrm{int},k}^{\prime}(x_{1}-x_{k})\bigg{)}\\\
&=\frac{1}{\gamma_{1}}\bigg{(}-\kappa_{\mathrm{ext}}(x_{1}-v_{0}t)+2\pi\sum_{k=2}^{d}u_{2k-2}u_{2k-3}\sin\big{(}2\pi
u_{2k-2}(x_{1}-x_{k})\big{)}\bigg{)}\\\
a_{i}(x,u,t)&=\frac{1}{\gamma_{i}}\,V_{\mathrm{int},i}^{\prime}(x_{1}-x_{i})=-\frac{2\pi}{\gamma_{i}}\,u_{2i-2}u_{2i-3}\sin\big{(}2\pi
u_{2i-2}(x_{1}-x_{i})\big{)}\end{aligned}$ for $i=2,\ldots,d$ and (106b)
$\displaystyle b_{i}(x,u,t)$ $\displaystyle=1/\gamma_{i},\qquad\qquad
i=1,\ldots,d.$

Notice that the $b_{i}$’s are constant and therefore independent of $(x,u,t)$.
Moreover, we have that

(107) $\displaystyle
a_{x}(x,u,t)=\left(\left(\partial_{x_{j}}a_{i}(x,u,t)\right)\right)\in{\mathbb{R}}^{d\times
d}$

with

$\displaystyle\begin{split}\partial_{x_{1}}a_{1}(x,u,t)&=-\frac{\kappa_{\mathrm{ext}}}{\gamma_{1}}+\frac{4\pi^{2}}{\gamma_{1}}\sum_{k=2}^{K+1}u_{2k-2}^{2}u_{2k-3}\cos\big{(}2\pi
u_{2k-2}(x_{1}-x_{k})\big{)},\\\
\partial_{x_{j}}a_{1}(x,u,t)&=-\frac{4\pi^{2}}{\gamma_{1}}\,u_{2k-2}^{2}u_{2k-3}\cos\big{(}2\pi
u_{2k-2}(x_{1}-x_{j})\big{)},\quad j=2,\ldots,K+1,\\\
\partial_{x_{1}}a_{i}(x,u,t)&=-\frac{4\pi^{2}}{\gamma_{i}}\,u_{2i-2}^{2}u_{2i-3}\cos\big{(}2\pi
u_{2i-2}(x_{1}-x_{i})\big{)},\quad i=2,\ldots,K+1,\\\
\partial_{x_{i}}a_{i}(x,u,t)&=\frac{4\pi^{2}}{\gamma_{i}}\,u_{2i-2}^{2}u_{2i-3}\cos\big{(}2\pi
u_{2i-2}(x_{1}-x_{i})\big{)},\quad i=2,\ldots,K+1.\end{split}$

We also find for

$\displaystyle
a_{u}(x,u,t)=\left(\left(\partial_{u_{j}}a_{i}(x,u,t)\right)\right)\in{\mathbb{R}}^{d\times
r},$

with the following components. For the right-hand side $a_{1}$ corresponding
to the tracer particle, we calculate

(108)
$\displaystyle\begin{split}\partial_{u_{2k-3}}a_{1}(x,u,t)&=\frac{2\pi}{\gamma_{1}}u_{2k-2}\sin\big{(}2\pi
u_{2k-2}(x_{1}-x_{k})\big{)},\\\
\partial_{u_{2k-2}}a_{1}(x,u,t)&=\frac{2\pi}{\gamma_{1}}\bigg{(}u_{2k-3}\sin\big{(}2\pi
u_{2k-2}(x_{1}-x_{k})\big{)}\\\ &\hskip 34.1433pt+2\pi
u_{2k-3}u_{2k-2}(x_{1}-x_{k})\cos\big{(}2\pi
u_{2k-2}(x_{1}-x_{k})\big{)}\bigg{)}.\end{split}$

For the other right-hand sides $a_{i}$, $i=2,\ldots,d$, we calculate

(109)
$\displaystyle\begin{split}\partial_{u_{2i-3}}a_{i}&=-\frac{2\pi}{\gamma_{i}}u_{2i-2}\sin\big{(}2\pi
u_{2i-2}(x_{1}-x_{i})\big{)},\\\
\partial_{u_{2i-2}}a_{i}&=-\frac{2\pi}{\gamma_{i}}\Big{(}u_{2i-3}\sin\big{(}2\pi
u_{2i-2}(x_{1}-x_{i})\big{)}\\\ &\qquad\qquad+2\pi
u_{2i-2}u_{2i-3}(x_{1}-x_{i})\sin\big{(}2\pi
u_{2i-2}(x_{1}-x_{i})\big{)}\Big{)}.\end{split}$

All other components vanish, i.e. $\partial_{u_{j}}a_{i}=0$ if $i\neq j$ for
$j=1,\ldots,d$, $i=2,\ldots,d$ since the position of each particle depends not
on parameters of the other particles. Hence, $\partial_{u}a(x,u,t)$ has the
structure

$\displaystyle\begin{pmatrix}\partial_{u_{1}}a_{0}&\partial_{u_{2}}a_{0}&\ldots&\ldots&\ldots&\ldots&\partial_{u_{2K}}a_{0}\\\
\partial_{u_{1}}a_{1}&0&0&\partial_{u_{K}}a_{1}&0&\ldots&0\\\
0&\ddots&0&0&\ddots&0&\vdots\\\ \vdots&0&\ddots&0&0&\ddots&0\\\
0&\ldots&0&\partial_{u_{K}}a_{K}&0&0&\partial_{u_{2K}}a_{N_{\mathrm{p}}}&\end{pmatrix}.$

Notice that $\partial_{u}b=0$ and $\partial_{x}b=0$. We define the following
continuous problem

(110a) $\displaystyle\min
J(\bm{X},u)=\frac{1}{2}\left\lVert[{\mathcal{C}}({\bm{X}})](\cdot)\,-\,{c^{\mathsf{d}}}(\cdot)\right\rVert_{L^{2}(0,T)}^{2}+\frac{1}{2}\left\lVert[{\mathcal{C}}({\bm{X}})](T)\,-\,{c^{\mathsf{d}}}(T)\right\rVert_{2}^{2},$
(110b) $\displaystyle\hskip 1.42262pt\text{s.t.
}(\bm{X},u)\in{\mathscr{X}}\times{\mathbb{R}}^{2{N_{\mathrm{p}}}}\text{ solves
\eqref{eq:Stochastic_Prantl_Tomlinson_intro}},$

where we have set

(111)
$\displaystyle[{\mathcal{C}}({\bm{X}})](t)\colonequals\frac{{\mathbb{E}}[Y_{1}(t)Y_{1}(0)]}{{\mathbb{E}}[Y_{1}^{2}(t)]}\quad\text{and}\quad
Y_{1}(t)=X_{1}-{\mathbb{E}}[X_{1}(t)],$

The term $[{\mathcal{C}}({\bm{X}})](t)$ stands for the (normalized)
correlation function and is connected to the _mean conditional displacement_
that is used in [25]. Notice that it holds for all $t$ that

$\displaystyle[{\mathcal{C}}(\bm{X})](t)=\frac{\mathrm{Cov}(Y_{1}(t),Y_{1}(0))}{{\mathbb{V}}[Y_{1}(t)]},$

since ${\mathbb{E}}[Y_{1}(t)]=0$ by the definition of $Y_{1}$.

######  Remark 7.4.

  1. 1)

The coefficients functions $a$ and $b$ are given by (106). Hence, they consist
of functions that are either affine linear in $x$ and $u$ or smooth and
bounded and therefore fulfill the Lipschitz property. Moreover, by linearity
and boundedness, they also fulfill (13). Consequently, Assumption 3.2 holds
also for the SPT model.

  2. 2)

Also the coefficients functions $a$ and $b$ are continuously differentiable
functions as a composition of smooth functions. Their first partial
derivatives are given by (108), (109), (7.2). Furthermore, the derivatives are
bounded in ${\bm{X}}$ due to the boundedness of the $\mathrm{sin}$ function
and bounded in $u$ due to our assumptions on ${{\mathscr{U}}_{\mathsf{ad}}}$.
Thus, Assumption 5.1-2) and Assumption 6.4 are fulfilled.

  3. 3)

We do not have the convexity of $\mathcal{C}$ since it contains the variance.
But assuming that the variance of $X_{u}$ is not vanishing, we have continuous
differentiability of the functional as a composition of continuous
differentiable functions and hence that Assumption 5.1-1) is fulfilled.
Furthermore, the functional is also bounded from below by zero, since it
consists of norms. Consequently, Assumption 4.1 is only partially satisfied.

Now, we formulate the discretized model and the discretized adjoint model to
the problem (110). For this, we define the discrete version of ${\mathcal{C}}$
as follows:

$\displaystyle{\mathcal{C}}_{\nu}^{M}({\bm{X}}^{h})\colonequals\frac{{\mathbb{E}}^{M}[Y_{1,\nu}^{h}Y_{1,0}^{h}]}{{\mathbb{E}}^{M}[(Y_{1,\nu}^{h})^{2}]}\quad\text{for
}\nu=1,\ldots,N\text{ and
}Y_{1,\nu}^{h}=X_{1,\nu}^{h}-{\mathbb{E}}^{M}\big{[}X_{1,\nu}^{h}\big{]}.$

Recall the definition of ${\mathbb{E}}^{M}$ in (51). Notice that the
derivative of the functional with respect to $X_{1,\nu}^{\mu}$,
$\nu\in\\{1,\ldots,N\\}$, $\mu\in\\{1,\ldots,M\\}$ is given by

$\displaystyle\left({\mathcal{C}}^{M}_{\nu}({\bm{X}}^{h})-{c_{\nu}^{\mathsf{d}}}\right)\bigg{(}\frac{1}{M^{2}}\sum_{\mu^{\prime}=1}^{M}(Y^{\mu^{\prime}}_{1,\nu})^{2}\,Y^{\mu}_{1,0}-\frac{1}{M^{2}}\sum_{\mu^{\prime}=1}^{M}\big{(}Y^{\mu^{\prime}}_{1,\nu}Y^{\mu^{\prime}}_{1,0}\big{)}\,2Y^{\mu}_{1,\nu}\bigg{)}\bigg{(}\frac{1}{M}\sum_{\mu^{\prime}=1}^{M}(Y^{\mu^{\prime}}_{1,\nu})^{2}\bigg{)}^{-2}$
$\displaystyle\quad=\frac{1}{M}\left({\mathcal{C}}^{M}_{\nu}({\bm{X}}^{h})-{c_{\nu}^{\mathsf{d}}}\right){\mathbb{E}}^{M}\big{[}(Y^{h}_{1,\nu})^{2}\big{]}\left(Y^{\mu}_{1,0}-2\,{\mathcal{C}}_{\nu}^{M}({\bm{X}}^{h})Y^{\mu}_{1,\nu}\right)\big{(}{\mathbb{E}}^{M}\big{[}(Y^{h}_{1,\nu})^{2}\big{]}\Big{)}^{-2}$
$\displaystyle\quad=\frac{1}{M}\left({\mathcal{C}}^{M}_{\nu}({\bm{X}}^{h})-{c_{\nu}^{\mathsf{d}}}\right)\frac{Y^{\mu}_{1,0}-2\,{\mathcal{C}}_{\nu}^{M}({\bm{X}}^{h})Y^{\mu}_{1,\nu}}{{\mathbb{E}}^{M}[(Y^{h}_{1,\nu})^{2}]}$

Thus, we can write the derivative of the functional with respect to the state
variable ${\bm{X}}^{h}=(X^{h}_{i})_{i=1,\ldots,d}$ as

(112) $\displaystyle
J^{h}_{{\bm{X}}^{h}}({\bm{X}}^{\mu}_{\nu})=\begin{cases}\frac{1}{M}\left({\mathcal{C}}^{M}_{\nu}({\bm{X}}^{h})-{c_{\nu}^{\mathsf{d}}}\right)\frac{Y^{\mu}_{1,0}-2\,{\mathcal{C}}_{\nu}^{M}({\bm{X}}^{h})Y^{\mu}_{1,\nu}}{{\mathbb{E}}^{M}[(Y^{h}_{1,\nu})^{2}]}&\text{
if }i=1,\\\ 0&\text{ otherwise}.\end{cases}$

We introduce the adjoint variable
${\bm{\Lambda}}^{h}=\\{{\bm{\Lambda}}^{\mu}_{\nu}\\}\in{\mathscr{X}}^{h}$.
Using (60a) and (107), we can formulate the discrete adjoint equation that
evolves backwards in time for $\nu=N-1,\ldots,0$

(113a) $\displaystyle{\bm{\Lambda}}^{\mu}_{\nu}=\left(I+\Delta
t\,a_{x}^{\nu}(\bm{X}^{\mu}_{\nu},u)\right)^{\top}{\bm{\Lambda}}^{\mu}_{\nu+1}+\Delta
t\,J^{h}_{{\bm{X}}^{h}}({\bm{X}}^{\mu}_{\nu})\quad\text{for }\mu=1,\ldots,m$
with the terminal condition (113b)
$\displaystyle{\bm{\Lambda}}^{\mu}_{N}=J^{h}_{{\bm{X}}^{h}}({\bm{X}}^{\mu}_{N})\quad\text{for
}\mu=1,\ldots,m.$

Furthermore, the discrete gradient is given by

(114)
$\displaystyle\nabla\hat{J}^{h}(\mathrm{u})=\frac{1}{M}\frac{1}{N-1}\sum_{\mu=1}^{M}\sum_{\nu=0}^{N-1}\left(a_{u}^{\nu}(\bm{X}^{\mu}_{\nu},u)\right)^{\top}{\bm{\Lambda}}_{\nu+1}^{\mu}.$

For the first numerical experiment using the SPT model, we consider a single
bath particle, i.e. we set $K=1$. In Fig. 2, we present the results. The plot
in Fig. 2(a) shows the simulation of (104) using the initial guess of the
parameters and the simulation using the parameters obtained by Algorithm 1.

(a) $K=1$

(b) $K=2$

Figure 2:  Results of numerical experiments for (110). Simulation of $X_{u}$
defined in (111) using the initial (grey) and optimal (black) guess of
parameters for (a) $K=1$ and (b) $K=2$, each time compared with the desired
${c^{\mathsf{d}}}$ (red).

For a second test, we consider (110) with two bath-particles, i.e. $K=2$. In
this case, we observed that a longer calibration time is needed, also because
the equilibration time $t_{eq}$ is longer (cf. Remark 7.3). We plot the
results in Fig. 2(b).

It is evident, that our method manages to find optimal parameters in order to
match the desired behavior of the model.

##  8 Acknowledgments

This work was partially funded by the Deutsche Forschungsgemeinschaft within
SFB1432, Project-ID 425217212. We would like to express our gratitude to
Juliana Caspers, Luis Reinalter as well as Clemens Bechinger and Matthias
Krüger for very helpful physical insights in the SPT model. We also thank
Behzad Azmi for very fruitful discussions.

## References

  * [1] Annunziato, M., and Borzì, A. A Fokker-Planck control framework for stochastic systems.  EMS Surv. Math. Sci. 5, 1-2 (2018), 65–98.
  * [2] Bartsch, J., and Borzì, A. MOCOKI: A Monte Carlo approach for optimal control in the force of a linear kinetic model.  Comput. Phys. Commun. 266 (2021), 108030.
  * [3] Bartsch, J., Nastasi, G., and Borzì, A. Optimal control of the Keilson-Storer master equation in a Monte Carlo framework.  J. Comput. Theor. Transp. 50, 5 (2021), 454–482.
  * [4] Bismut, J.-M. Conjugate convex functions in optimal stochastic control.  J. Math. Anal. Appl. 44 (1973), 384–404.
  * [5] Bismut, J.-M. An introductory approach to duality in optimal stochastic control.  SIAM Rev. 20, 1 (1978), 62–78.
  * [6] Bottou, L., Curtis, F. E., and Nocedal, J. Optimization methods for large-scale machine learning.  SIAM Rev. 60, 2 (2018), 223–311.
  * [7] Breitenbach, T., and Borzí, A. The Pontryagin maximum principle for solving Fokker-Planck optimal control problems.  Comput. Optim. Appl. 76, 2 (2020), 499–533.
  * [8] Buckdahn, R., Labed, B., Rainer, C., and Tamer, L. Existence of an optimal control for stochastic control systems with nonlinear cost functional.  Stochastics 82, 3 (2010), 241–256.
  * [9] Clevenhaus, A., Totzeck, C., and Ehrhardt, M. A gradient based calibration method for the heston model.  arXiv (2023).
  * [10] Conlan, J., and Wang, C.-L. A unified approach to continuous and discrete Gronwall-Bellman inequalities.  Appl. Anal. 44, 3-4 (1992), 243–252.
  * [11] Cresson, J., and Sonner, S. A note on a derivation method for SDE models: applications in biology and viability criteria.  Stoch. Anal. Appl. 36, 2 (2018), 224–239.
  * [12] El Karoui, N., Hu̇u̇ Nguyen, D., and Jeanblanc-Picqué, M. Compactification methods in the control of degenerate diffusions: existence of an optimal control.  Stochastics 20, 3 (1987), 169–219.
  * [13] Estevao, V. M., and SaÈrndal, C.-E. A functional form approach to calibration.  J Off Stat 16, 4 (2000), 379.
  * [14] Evans, L. C. An introduction to stochastic differential equations.  American Mathematical Society, Providence, RI, 2013.
  * [15] Fabbri, G., Gozzi, F., and Świećh, A. Stochastic optimal control in infinite dimension, vol. 82 of Probability Theory and Stochastic Modelling.  Springer, Cham, 2017.  Dynamic programming and HJB equations, With a contribution by Marco Fuhrman and Gianmario Tessitore.
  * [16] Fabrizio, E., and Monetti, V. Methodologies and advancements in the calibration of building energy models.  Energies 8, 4 (2015), 2548–2574.
  * [17] Gilli, M., and Schumann, E. Calibrating option pricing models with heuristics.  In Natural computing in computational finance. Springer, 2011, pp. 9–37.
  * [18] Gong, B., Liu, W., Tang, T., Zhao, W., and Zhou, T. An efficient gradient projection method for stochastic optimal control problems.  SIAM J Numer Anal 55, 6 (2017), 2982–3005.
  * [19] Graham, C., and Talay, D. Stochastic simulation and Monte Carlo methods, vol. 68 of Stochastic Modelling and Applied Probability.  Springer, Heidelberg, 2013.  Mathematical foundations of stochastic simulation.
  * [20] Groß, B. Applications of the Adjoint Method in Stochastic Financial Modelling.  PhD thesis, University of Trier, 2015.
  * [21] Higham, D. J. An algorithmic introduction to numerical simulation of stochastic differential equations.  SIAM Rev. 43, 3 (2001), 525–546.
  * [22] Hinze, M., Pinnau, R., Ulbrich, M., and Ulbrich, S. Optimization with PDE constraints, vol. 23 of Mathematical Modelling: Theory and Applications.  Springer, New York, 2009.
  * [23] Jagla, E. A. The Prandtl-Tomlinson model of friction with stochastic driving.  J. Stat. Mech. Theory Exp., 1 (2018), 013401, 14.
  * [24] Jain, R., Ginot, F., Berner, J., Bechinger, C., and Krüger, M. Two step micro-rheological behavior in a viscoelastic fluid.  J. Chem. Phys. 154, 18 (2021), 184904.
  * [25] Jain, R., Ginot, F., and Krüger, M. Micro-rheology of a particle in a nonlinear bath: Stochastic Prandtl-Tomlinson model.  Phys. Fluids 33, 10 (Oct. 2021), 103101.
  * [26] Kaebe, C., Maruhn, J. H., and Sachs, E. W. Adjoint-based Monte Carlo calibration of financial market models.  Finance Stoch. 13, 3 (2009), 351–379.
  * [27] Karatzas, I., and Shreve, S. E. Brownian motion and stochastic calculus, vol. 113 of Graduate Texts in Mathematics.  Springer-Verlag, New York, 1988.
  * [28] Kawasaki, K. Simple derivations of generalized linear and nonlinear Langevin equations.  J. Phys. A Math. Theor. 6, 9 (Sept. 1973), 1289–1295.
  * [29] Kosmol, P., and Pavon, M. Lagrange approach to the optimal control of diffusions.  Acta Appl. Math. 32, 2 (1993), 101–122.
  * [30] Kubo, R. The fluctuation-dissipation theorem.  Rep. Prog. Phys. 29, 1 (1966), 255.
  * [31] Ló, B. P., Haslam, A., and Adjiman, C. An algorithm for the estimation of parameters in models with stochastic differential equations.  Chem. Eng. Sci. 63, 19 (2008), 4820–4833.
  * [32] Lü, Q., and Zhang, X. General Pontryagin-type stochastic maximum principle and backward stochastic evolution equations in infinite dimensions.  SpringerBriefs in Mathematics. Springer, Cham, 2014.
  * [33] Luenberger, D. G. Optimization by Vector Space Methods.  John Wiley & Sons, Inc., New York, NY, 1969.
  * [34] Ma, J., and Yong, J. Forward-backward stochastic differential equations and their applications, vol. 1702 of Lecture Notes in Mathematics.  Springer-Verlag, Berlin, 1999.
  * [35] Mao, X. Stochastic differential equations and applications, second ed.  Horwood Publishing Limited, Chichester, 2008.
  * [36] Maruyama, G. Continuous Markov processes and stochastic equations.  Rend. Circ. Mat. Palermo (2) 4 (1955), 48–90.
  * [37] Milstein, G. N., and Tretyakov, M. V. Stochastic numerics for mathematical physics, second ed.  Scientific Computation. Springer, Cham, [2021] ©2021.
  * [38] Müller, B., Berner, J., Bechinger, C., and Krüger, M. Properties of a nonlinear bath: experiments, theory, and a stochastic prandtl–tomlinson model.  New J. Phys. 22, 2 (2020), 023014.
  * [39] Øksendal, B. Stochastic differential equations: an introduction with applications.  Springer Science & Business Media, 2013.
  * [40] Pham, H. Continuous-time stochastic control and optimization with financial applications, vol. 61 of Stochastic Modelling and Applied Probability.  Springer-Verlag, Berlin, 2009.
  * [41] Protter, P. E. Stochastic integration and differential equations, second ed., vol. 21 of Stochastic Modelling and Applied Probability.  Springer-Verlag, Berlin, 2005.  Corrected third printing.
  * [42] Renardy, M., and Rogers, R. C. An introduction to partial differential equations, second ed., vol. 13 of Texts in Applied Mathematics.  Springer-Verlag, New York, 2004.
  * [43] Sachs, E. W., and Schu, M. Gradient computation for model calibration with pointwise observations.  In Control and optimization with PDE constraints, vol. 164 of Internat. Ser. Numer. Math. Birkhäuser/Springer Basel AG, Basel, 2013, pp. 117–136.
  * [44] Tröltzsch, F. Optimal control of partial differential equations: Theory, methods, and applications, vol. 112 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, Providence, RI, 2010.
  * [45] van Kampen, N. G. Stochastic differential equations.  Phys. Rep. 24, 3 (1976), 171–228.
  * [46] van Kampen, N. G. Stochastic processes in physics and chemistry, vol. 888 of Lecture Notes in Mathematics.  North-Holland Publishing Co., Amsterdam-New York, 1981.
  * [47] Yong, J., and Zhou, X. Y. Stochastic Controls, vol. 43 of Applications of Mathematics (New York).  Springer-Verlag, New York, 1999.  Hamiltonian systems and HJB equations.
  * [48] Zwanzig, R. Nonlinear generalized Langevin equations.  J. Stat. Phys. 9, 3 (1973), 215–220.