#  Anatomy of linear and non-linear intermolecular exchange in $S=1$
nanographenes

J. C. G. Henriques1,2, J. Fernández-Rossier1111On permanent leave from
Departamento de Física Aplicada, Universidad de Alicante, 03690 San Vicente
del Raspeig, Spain. (1) International Iberian Nanotechnology Laboratory (INL),
Av. Mestre José Veiga, 4715-330 Braga, Portugal  (2) Universidade de Santiago
de Compostela, 15782 Santiago de Compostela, Spain

###### Abstract

Nanographene triangulenes with a $S=1$ ground state have been used as building
blocks of antiferromagnetic Haldane spin chains realizing a symmetry protected
topological phase. By means of inelastic electron spectroscopy, it was found
that the intermolecular exchange contains both linear and non-linear
interactions, realizing the bilinear-biquadratic Hamiltonian. Starting from a
Hubbard model, and mapping it to an interacting Creutz ladder, we analytically
derive these effective spin-interactions using perturbation theory, up to
fourth order. We find that for chains with more than two units other
interactions arise, with same order-of-magnitude strength, that entail second
neighbor linear, and three-site non-linear exchange. Our analytical
expressions compare well with experimental and numerical results. We discuss
the extension to general $S=1$ molecules, and give numerical results for the
strength of the non-linear exchange for several nanographenes. Our results
pave the way towards rational design of spin Hamiltonians for nanographene
based spin chains.

Non-linear exchange, i.e., spin interactions that go beyond the simple
Heisenberg coupling between two spins, play a prominent role in many physical
systems, such as antiferromagnetic transition metal oxides[1], magnetic
impurities in insulators[2], magnetic multilayers[3], chiral magnets[4] and
magnetic two-dimensional materials[5]. Non-linear exchange is a key ingredient
in the exactly solvable AKLT models[6], whose ground state is a resource for
measurement based quantum computing[7].

The relative size and sign of linear and non-linear exchange can have a
dramatic impact in several contexts. For instance, whereas the swap gate, or
permutation operator, for spin qubits can be implemented with a linear
Heisenberg interaction[8], for spin qudits it requires the presence of non-
linear exchange terms[9]. Alternatively, in the case of the two-dimensional
$S=3/2$ honeycomb lattice, the relative size of linear and non-linear exchange
controls the nature of the ground state and its excitation spectrum [6, 10,
11].

Here we undertake the exploration of non-linear exchange in $S=1$
nanographenes. This class of system features outstanding flexibility to
realize $S=1$ molecules with different shapes and sizes. It has been recently
shown[12] that there are 383 different nanographenes that can be formed with 9
hexagons or less. Therefore, this type of system provides an ideal arena to
engineer intermolecular exchange. One of the simplest $S=1$ nanographenes is
the so-called [3]-triangulene.

Triangulenes are graphene fragments with the shape of an equilateral triangle,
of various sizes and terminated with zigzag edges; these are customarily
defined in terms of the number of benzenes, $n$, in a given edge[13, 14]
(termed a [n]-triangulene). Based on Lieb’s theorem for the Hubbard model for
bipartite lattices at half-filling [15], $[n]$-triangulenes are predicted to
be open-shell multiradicals, with the spin of the ground state scaling as
$2S=n-1$, associated with a half-full shell of $n-1$ in-gap non-bonding zero
modes[16, 13, 17, 18, 19, 20].

Due to recent breakthroughs in bottom-up synthesis techniques [21, 22, 23, 24,
25], and the capability of atomic precision manipulation of organic molecules,
triangulenes have been used as building blocks of larger molecular structures
[26, 27, 28, 29]. A prime example of this is the recent realization of Haldane
spin chains, where [3]-triangulenes (henceforth referred simply as
triangulene) were coupled in order to generate chains with more than 16 units
[27]. There, the inelastic electron tunneling spectroscopy (IETS) was
described with an effective $S=1$ spin Hamiltonian that included both linear
and non-linear exchange terms, i.e. the so called BLBQ Hamiltonian [30, 31,
32, 33, 34]:

$\displaystyle
H_{\text{BLBQ}}=\sum_{i}J\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{i+1}+B\left(\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{i+1}\right)^{2},$
(1)

where each triangulene is represented by a spin-1 operator
$\boldsymbol{S}_{i}$, the sum runs over all triangulenes in the chain and the
parameters $J$ and $B$ are the linear and non-linear exchange couplings,
respectively. The introduction of the non-linear exchange term proved
essential to increase the accuracy of the spin model when compared with
experimental data [27] and full fermionic numerical approaches [32].

In the following, we derive the effective $S=1$ spin Hamiltonian (1) starting
from a fermionic model for the nanographenes. Importantly, our derivation
unveils the presence of additional second neighbor linear interactions, and
non-linear exchange that involve three-spin terms, with strength comparable to
$B$ in Eq. (1). As a fermionic model we use a single-orbital Hubbard model
(see Supporting Information), containing hoppings between first ($t$) and
third neighbor ($t_{3}$) sites, and an on-site Hubbard repulsion term ($U$)
which deals with the intra-atomic Coulomb repulsion cost associated with
having a given $\pi$-orbital doubly occupied [20]. The Hubbard model has been
validated by comparison with multiconfigurational calculations obtained with
full-quantum chemistry _ab-initio_ methods [20], for $t=-2.7$ eV, $t_{3}\sim
t/10$ and $U\sim|t|$ [35].

As a starting point, we employ a Hubbard Hamiltonian at half-filling to
describe the smallest chain, i.e. a triangulene dimer. In order to derive the
effective spin model, we describe the fermionic many-body states of the
nanographenes in a truncated Hilbert space, corresponding to the complete
active space (CAS) approximation (see Supporting Information) [20, 27, 36,
37]. In Fig. 1 (a) and (b) we show the single-particle energies for the
triangulene dimer with $t_{3}=t/10$ and $t_{3}=0$; in both cases two possible
choices of active spaces are indicated. When $t_{3}=0$ one finds four states
with zero energy; these correspond to the unhybridized zero modes of the
individual triangulenes. For $t_{3}=t/10$, intermolecular hybridization
promoted by $t_{3}$ lifts this degeneracy. In panels (c) and (d) of the same
figure, the absolute value of the site representation of the four modes
closest to zero energy are depicted for $t_{3}=t/10$ and $t_{3}=0$. For the
latter case, the wave functions were chosen as eigenfunctions of the $C_{3}$
symmetry operator with eigenvalues $\pm 2\pi/3$. For future reference, we
refer to panel (c) as the molecular orbital basis, and to panel (d) as the
$C_{3}$ symmetric basis.

Figure 1: Single particle energies for the triangulene dimer with $t_{3}=t/10$
(a) and $t_{3}=0$ (b). Energies are given in unit of $|t|=2.7$ eV. The dashed
boxes indicate two possible choices for the active space where the Hubbard
Hamiltonian is represented. CAS($N_{O},N_{e}$) refers to an active space with
$N_{O}$ modes where $N_{e}$ electrons are distributed). Panels (c) and (d)
show the absolute value of the site representation of the four molecular
orbitals closer to zero energy for $t_{3}=t/10$ and $t_{3}=0$, respectively.
Each circle is colored according to the sublattice where it is located.

To gain physical insight about the properties of this system, we represent the
Hamiltonian in the $C_{3}$ symmetric basis taking into account only the four
modes at zero energy. The representation of the Hubbard model in that basis
leads to the following effective Hamiltonian:

$\displaystyle H$
$\displaystyle=\sum_{\mu^{\prime}\mu}\sum_{\sigma}\tau_{\mu^{\prime}\mu}d^{\dagger}_{\mu^{\prime},\sigma}d_{\mu,\sigma}+\sum_{\mu}\mathcal{U_{\mu}}n_{\mu,\uparrow}n_{\mu,\downarrow}$
$\displaystyle+\sum_{\triangle}\mathcal{J}_{\triangle_{+},\triangle_{-}}\left(n_{\triangle_{+},\uparrow}n_{\triangle_{-},\downarrow}+n_{\triangle_{-},\uparrow}n_{\triangle_{+},\downarrow}\right)$
$\displaystyle+\sum_{\triangle}\mathcal{J}_{\triangle_{+},\triangle_{-}}\left(d^{\dagger}_{\triangle_{+},\uparrow}d_{\triangle_{-},\uparrow}d^{\dagger}_{\triangle_{-},\downarrow}d_{\triangle_{+},\downarrow}+\text{h.c.}\right)$
(2)

where $\mu$ and $\mu^{\prime}$ run over the four $C_{3}$ symmetric modes and
$\sigma=\uparrow,\downarrow$. The sums over $\triangle$ run over the two
triangulenes, and $\triangle_{\pm}$ refers to the modes with eigenvalues $\pm
2\pi/3$ in a given triangulene. The operator $d^{(\dagger)}_{\mu,\sigma}$
annihilates (creates) an electron in the $\mu$ mode with spin $\sigma$, and
$n_{\mu,\sigma}=d^{\dagger}_{\mu,\sigma}d_{\mu,\sigma}$.

This Hamiltonian contains three distinct types of terms. The first term in Eq.
(2) represents the hopping between two modes, $\mu$ and $\mu^{\prime}$, with
an amplitude $\tau_{\mu\mu^{\prime}}=t_{3}\sum_{\langle\langle\langle
i,j\rangle\rangle\rangle}\Phi^{*}_{\mu}(i)\Phi_{\mu^{\prime}}(j)$ where
$\Phi_{\mu}(i)$ is the site representation of the $\mu$-mode (depicted in Fig.
1). This hopping amplitude is zero when the two modes are localized in the
same sublattice; thus, hoppings between modes in the same triangulene vanish.
In addition, due to the $C_{3}$ symmetry of the $\mu$-modes, all the finite
hoppings have the same absolute value, which we label simply as $\tau$. The
pictorial representation of this hopping structure, shown in fig. 2, reveals
the single-particle Hamiltonian maps into a Creutz ladder model[38] with
vanishing vertical hoppings.

The remaining terms in the Hamiltonian account for electron-electron
interactions. The second term of Eq. (2), corresponds to an effective Hubbard
repulsion that penalizes the double occupancy of a given mode $\mu$, with the
energy penalty being proportional to the inverse participation ratio (IPR),
and is defined as $\mathcal{U}_{\mu}=U\sum_{i}|\Phi_{\mu}(i)|^{4}$, where the
sum runs over all the sites. Again, due to the $C_{3}$ symmetry of the modes,
one finds $\mathcal{U}_{\mu}\equiv\mathcal{U}$ is independent of $\mu$.

The last two lines in Eq. (2) describe an intra-triangulene ferromagnetic
exchange interaction, with the exchange coupling given by
$\mathcal{J}_{\mu\mu^{\prime}}=U\sum_{i}|\Phi_{\mu}(i)|^{2}|\Phi_{\mu^{\prime}}(i)|^{2}$.
As before, $\mathcal{J}_{\mu\mu^{\prime}}=\mathcal{J}$ is independent of $\mu$
and $\mu^{\prime}$, and $\mathcal{U}=\mathcal{J}$.

Figure 2: (a) Schematic representation of the interactions present in the
Hamiltonian of Eq. (2), with modes in the same triangulene ferromagnetically
coupled, and hoppings between the two triangulenes; (b) Extension of the model
for the case of a triangulene chain.

Therefore the model (2) realizes a Creutz-ladder [38], with the difference
that the vertical hoppings are replaced by a ferromagnetic exchange and there
are Hubbard interactions. In Fig. 2 we present a pictorial representation of
the Hamiltonian of Eq. (2) for the triangulene dimer, as well as the extension
to larger structures, e.g. a triangulene chain. The extension of Eq. (2) to
triangulene chains is easily obtained by running the sums over all
triangulenes (see Supporting Information for the numerical solution of Eq.(2)
for a triangulene dimer and trimer).

Let us note in passing, that for molecules less symmetric than triangulenes,
for example panels (c) to (f) of Fig. 4, the hoppings
$\tau_{\mu\mu^{\prime}}$, the effective Hubbard repulsion $\mathcal{U}_{\mu}$,
and the exchange coupling $\mathcal{J}_{\mu\mu^{\prime}}$ are in principle
mode-dependent. Moreover, it might be necessary to include a new term
associated with electron-pair hoppings in Eq. (2). As discussed in [20], this
term vanishes for triangulenes, but may be finite for other molecules.
Including it is straightforward, and does not affect the main physical
features we presently wish to discuss.

Starting from the many-body fermionic model just presented, we now produce an
effective low energy description which can be related to the effective BLBQ
spin model. To this end, we shall employ degenerate perturbation up to fourth
order, treating the hoppings $\tau$ as a perturbation.

For the unperturbed system ($\tau=0$) the two pairs of modes are decoupled,
and each triangulene can be studied individually. At half-filling, with two
electrons per triangulene, the lowest energy states in a given triangulene are
$|\uparrow\rangle|\uparrow\rangle$, $|\downarrow\rangle|\downarrow\rangle$ and
$\left(|\uparrow\rangle|\downarrow\rangle+|\downarrow\rangle|\uparrow\rangle\right)/\sqrt{2}$,
where each ket refers to one of the $C_{3}$ symmetric modes in a triangulene.
These states correspond to the three spin projections of a $S=1$ state formed
by the ferromagnetic coupling of two spin-1/2 electrons, i.e.
$|\uparrow\rangle|\uparrow\rangle\equiv|\Uparrow\rangle$,
$|\downarrow\rangle|\downarrow\rangle\equiv|\Downarrow\rangle$ and
$\left(|\uparrow\rangle|\downarrow\rangle+|\downarrow\rangle|\uparrow\rangle\right)/\sqrt{2}\equiv|+\rangle$,
with the ket on the right hand side referring to the state of the triangulene
as a whole. With these three states in each triangulene, the total unperturbed
dimer Hamiltonian at half-filling is nine-fold degenerate.

Representing Eq. (1) in the basis of two spin-1 objects, and focusing on the
$S_{z}=0$ sector, one finds that while the bilinear term, proportional to $J$,
is responsible for connecting the state $|+\rangle|+\rangle$ with the states
$|\Uparrow\rangle|\Downarrow\rangle$ and $|\Downarrow\rangle|\Uparrow\rangle$,
the biquadratic term, proportional to $B$, unlocks a new interaction between
the states $|\Uparrow\rangle|\Downarrow\rangle$ and
$|\Downarrow\rangle|\Uparrow\rangle$. Recalling the definition of the spin-1
states in terms of two spin-1/2, one realizes that while the processes
mediated by $J$ link states which differ by two spin flips, the process
mediated by $B$ connects states differing by four spin flips. Hence, the
biquadratic interaction is only to be expected in 4th order perturbation
theory, while the bilinear term should already be present in 2nd order.

The expressions for the 2nd and 4th order corrections in degenerate
perturbation theory read [39, 40]:

$\displaystyle h^{(2)}_{kp}$ $\displaystyle=\sum_{m}\frac{\langle
k|H_{\tau}|m\rangle\langle m|H_{\tau}|p\rangle}{\Delta E^{(0)}_{m}}$ (3)
$\displaystyle h^{(4)}_{kp}$
$\displaystyle=\sum_{mm^{\prime}m^{\prime\prime}}\frac{\langle
k|H_{\tau}|m\rangle\langle m|H_{\tau}|m^{\prime}\rangle\langle
m^{\prime}|H_{\tau}|m^{\prime\prime}\rangle\langle
m^{\prime\prime}|H_{\tau}|p\rangle}{\Delta E^{(0)}_{m}\Delta
E^{(0)}_{m^{\prime}}\Delta E^{(0)}_{m^{\prime\prime}}}$
$\displaystyle-\sum_{mm^{\prime}l}\frac{\langle k|H_{\tau}|m\rangle\langle
m|H_{\tau}|l\rangle\langle l|H_{\tau}|m^{\prime}\rangle\langle
m^{\prime}|H_{\tau}|p\rangle}{\Delta E^{(0)}_{m}\Delta E^{(0)}_{m}\Delta
E^{(0)}_{m^{\prime}}}$ (4)

where $k$, $p$ and $l$ label states inside the degenerate ground state, and
the sums run over all the states outside that subspace. The perturbation is
$H_{\tau}=\sum_{\mu\mu^{\prime}}\sum_{\sigma}\tau_{\mu\mu^{\prime}}d^{\dagger}_{\mu\sigma}d_{\mu\sigma}$
and $\Delta E^{(0)}_{m}$ is the unperturbed energy of the $|m\rangle$ state
relative to the ground state. In the absence of applied magnetic field, the
odd-order corrections vanish identically [41]. One important aspect to note is
that even though the initial and final state correspond to open-shell
configurations with two electrons per triangulene, the intermediate states
contain closed-shell configurations, as well as charge excitations, with
different number of electrons in the two triangulenes.

In 2nd order perturbation theory, where only two electron flips are
considered, a finite contribution for $J$ is found, which for the present
system simply reads $J^{(2)}=2\tau^{2}/\mathcal{U}$ (see Supporting
Information). This result is similar to the one usually found when when a
Hubbard chain is mapped to a Heisenberg chain of antiferromagnetically coupled
spins; the different numerical pre-factor steming from the different geometry
of our system.

Progressing to 4th order perturbation theory, where processes involving up to
four electron are considered, different paths appear which allow for
$|\Uparrow\rangle|\Downarrow\rangle\rightarrow|\Downarrow\rangle|\Uparrow\rangle$.
Carrying out the necessary calculations, one finds that finite contributions
appear for both $J$ and $B$. These are given by
$J^{(4)}=4\tau^{4}/\mathcal{U}^{3}$ and $B^{(4)}=8\tau^{4}/\mathcal{U}^{3}$.
The 2nd order contribution to $J$ dominates its 4th order counterpart in the
physically relevant region of the parameter space, and $J^{(4)}$ may be
neglected. Hence, combining the results from second and forth order
perturbation theory, we find in leading order:

$\displaystyle J\approx\frac{2\tau^{2}}{\mathcal{U}},\quad
B\approx\frac{8\tau^{4}}{\mathcal{U}^{3}},\quad\beta\equiv\frac{B}{J}\approx\frac{4\tau^{2}}{\mathcal{U}^{2}}.$
(5)

Using the definitions of $\tau$ [35] and $\mathcal{U}$ [20], and assuming
$t_{3}\approx t/10$, one can estimate the strength of the quadratic exchange
with respect to the linear one, characterized by $\beta$, using $\beta\approx
0.1(t/U)^{2}$. This rough approximation yields a nice agreement with the value
$\beta=0.09$ that was found in the description of experimental data in [27].
This good agreement, together with the analytical expressions of Eq. (5),
indicate that non-linear exchange is a higher-order manifestation of the same
underlying kinetic exchange mechanism[1] that gives rise to linear exchange.

To validate our analytical expressions we shall now compare them with the
results found from numerical diagonalization. The numerical results are
obtained by first diagonalizing Eq. (2), followed by matching the energies of
the first excitations with those of the BLBQ Hamiltonian (see Supporting
Information for details). This gives numerically both $J$ and $B$ as a
function of the parameters of the microscopic model, $U$ and $t_{3}$, and
allows the comparison with Eq. (5).

In Fig. 3 we show the comparison between the analytical expressions found for
$J$ and $B$ and the numerical results.

Figure 3: Comparison between the analytical expression and the result found
from numerical diagonalization for (a) the linear exchange $J$ and the non-
linear exchange $B$, and (b) the ratio $\beta=B/J$ as a function of the on-
site Hubbard repulsion $U$ in units of $|t|=2.7$ eV with $t_{3}=t/10$.

Inspecting this figure reveals an excellent agreement between the two
approaches, especially in the region $U\gg|t|$ where perturbation theory is
valid. Crucially, this agreement holds near $U\sim|t|$, the physically
relevant region. For $U=|t|$, one finds the linear exchange $J=20$ meV and
$\beta=0.11$. As $U$ increases $\beta$ approaches 0, and a Heisenberg-like
picture is recovered. For $U<<|t|$, outside the validity region of
perturbation theory, the mapping to a spin model fails since the order of low
energy excitations is no longer the same in the two approaches. The fact that
$\beta$ is bounded by the AKLT limit $\beta=1/3$ is a consequence of Lieb’s
theorem, which prevents a four-fold degenerate ground state.

Having completed the study of the triangulene dimer, we extend our analysis to
larger chains. Generalizing Eq. (2) to include more triangulenes, and once
again employing perturbation theory up to 4th order, we find that for a chain
composed of $N$ triangulenes the effective spin Hamiltonian reads (see
Supporting Information for detailed derivation):

$\displaystyle H_{N}$
$\displaystyle=J\sum_{i=1}^{N-1}\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{i+1}+B\sum_{i=1}^{N-1}(\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{i+1})^{2}$
$\displaystyle+J_{2}\sum_{i=1}^{N-2}\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{i+2}$
$\displaystyle+B_{1,1}\sum_{i=1}^{N-2}(\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{i+1})(\boldsymbol{S}_{i+1}\cdot\boldsymbol{S}_{i+2})+\text{h.c.}.$
(6)

While the first line corresponds to the BLBQ model for $N$ triangulenes, the
second and third lines contain new exchange interactions; the former describes
an antiferromagnetic second neighbor linear exchange, and the latter a
ferromagnetic quadratic exchange involving two adjacent triangulene pairs. We
note that the antiferromagnetic second neighbor linear exchange might promote
frustration if its strength becomes comparable with $J$. If perturbation
theory is extended up to sixth or eight order, additional exchange
interactions appear; however, since these are in higher order of
$\tau/\mathcal{U}$ we neglect their contribution. Crucially, the terms
proportional to $J_{2}$ and $B_{1,1}$ appear in the same order of perturbation
theory as $B$ and therefore must be accounted for. To leading order, we find:

$\displaystyle
J=\frac{2\tau^{2}}{\mathcal{U}},B=\frac{8\tau^{4}}{\mathcal{U}^{3}},J_{2}=\frac{79\tau^{4}}{12\mathcal{U}^{3}},B_{1,1}=-\frac{37\tau^{4}}{12\mathcal{U}^{3}}$
(7)

Hence, a consistent description of triangulene spin chains should include this
terms, missing in previous analysis[27].

We now briefly explore the non-linear exchange of various $S=1$ dimers. By
applying numerical diagonalization in the minimal CAS, and comparing it with
the excitation energies of the BLBQ Hamiltonian, we obtain the value of
$\beta$ for different dimers composed of molecules with a triplet ground
state. This should predict the relevance of non-linear exchange, as well as
its tunability, in chains formed with different building blocks. In Fig. 4 (b)
we show a triangulene dimer, where a benzene is introduced in the middle.
Formally, this molecule is nearly identical to the dimer we considered in the
text (depicted in panel (a)), with the main difference being that due to the
extra benzene, intermolecular hybridization between the two triangulenes is
greatly diminished, resulting in a smaller $\tau$. This leads to
$\beta=0.002$, two orders of magnitude smaller than what we found for the case
without the spacer (see Supporting Information for details on the change in
this value due to a larger active space). Next, we consider the molecule of
panel (c), where each monomer is composed of two Phenalenyl side-by-side; for
this dimer $\beta=0.14$. As benzenes are added between the two Phenalenyl,
depicted in panels (d) [42] and (e), the value of $\beta$ increases to
$\beta=0.16$ and $\beta=0.18$, respectively. This increase in $\beta$ may be
ascribed to the decrease of the IPR of the zero modes of the individual
molecules. At last, we consider panel (f), where an additional benzene is
added close to the binding site of the two molecules. This leads to a
significant decrease of $\beta$ to $\beta=0.03$. The reason for this sharp
decrease is similar to the one we gave when discussing panel (b). Hence, we
see that by carefully choosing the geometry of the molecules used as building
blocks, it should be possible to engineer the strength of non-linear exchange.

Figure 4: Values obtained for $\beta=B/J$ by numerically diagonalizing the
Hubbard Hamiltonian in a restricted Hilbert space with 4 modes only, for
dimers obtained from $S=1$ building blocks, taking $U=|t|$, $t_{3}=t/10$ and
$|t|=2.7$eV.

Even though we restricted our analysis to the minimal CAS, accounting only for
the zero modes, additional orbitals could have been included in the
calculation (see Fig. 1). This would introduce the so-called Coulomb driven
super-exchange discussed in [37]. We have verified that including this
additional mechanism would slightly change the numerical results, but preserve
the qualitative features of the model.

_Conclusions:_ We have considered chains of $S=1$ nanographenes, taking the
case of triangulene chains as the prototypical example. Using perturbation
theory we have found that non-linear exchange interactions are a higher-order
manifestation of the same mechanisms that give rise to linear exchange, and we
have obtained analytical expressions for their amplitude. The analytical
results permit us to relate molecule geometry with $\beta$, the degree of
exchange non-linearity. In addition, our analysis with more than two molecules
shows that new terms appear in the Hamiltonian, namely a second neighbor
linear and a three-site non linear exchange interaction, going beyond the BLBQ
paradigm. Future work will address the impact of these extra terms on the well
established phase diagram of the BLBQ model.

We acknowledge fruitful discussions with Gonçalo Catarina, António Costa and
David Jacob. We acknowledge financial support from FCT (Grant No. PTDC/FIS-
MAC/2045/2021), SNF Sinergia (Grant Pimag), FEDER /Junta de Andalucía, (Grant
No. P18-FR-4834), Generalitat Valenciana funding Prometeo2021/017 and
MFA/2022/045, and funding from MICIIN-Spain (Grant No. PID2019-109539GB-C41).

## References

  * Anderson [1959] P. W. Anderson, Physical Review 115, 2 (1959).
  * Harris and Owen [1963] E. Harris and J. Owen, Physical review letters 11, 9 (1963).
  * Slonczewski [1991] J. Slonczewski, Physical Review Letters 67, 3172 (1991).
  * Paul  _et al._ [2020] S. Paul, S. Haldar, S. Von Malottki, and S. Heinze, Nature communications 11, 4756 (2020).
  * Kartsev  _et al._ [2020] A. Kartsev, M. Augustin, R. F. Evans, K. S. Novoselov, and E. J. Santos, npj Computational Materials 6, 150 (2020).
  * Affleck  _et al._ [1987] I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. Rev. Lett. 59, 799 (1987).
  * Wei  _et al._ [2011] T.-C. Wei, I. Affleck, and R. Raussendorf, Phys. Rev. Lett. 106, 070501 (2011).
  * DiVincenzo  _et al._ [2000] D. P. DiVincenzo, D. Bacon, J. Kempe, G. Burkard, and K. B. Whaley, nature 408, 339 (2000).
  * Segraves [1964] P. H. Segraves,  _Representation of Permutation Operator in Quantum Mechanics_ , Master’s thesis, University of British Columbia (1964).
  * Pomata and Wei [2020] N. Pomata and T.-C. Wei, Phys. Rev. Lett. 124, 177203 (2020).
  * Ganesh  _et al._ [2011] R. Ganesh, D. Sheng, Y.-J. Kim, and A. Paramekanti, Physical Review B 83, 144414 (2011).
  * Yan  _et al._ [2023] Y. Yan, F. Zheng, B. Qie, J. Lu, H. Jiang, Z. Zhu, and Q. Sun, The Journal of Physical Chemistry Letters 14, 3193 (2023).
  * Fernández-Rossier and Palacios [2007] J. Fernández-Rossier and J. J. Palacios, Physical Review Letters 99, 177204 (2007).
  * Su  _et al._ [2020a] J. Su, M. Telychko, S. Song, and J. Lu, Angewandte Chemie International Edition 59, 7658 (2020a).
  * Lieb [1989] E. H. Lieb, Physical review letters 62, 1201 (1989).
  * Ovchinnikov [1978] A. A. Ovchinnikov, Theoretica Chimica Acta 47, 297 (1978).
  * Wang  _et al._ [2008] W. L. Wang, S. Meng, and E. Kaxiras, Nano letters 8, 241 (2008).
  * Wang  _et al._ [2009] W. L. Wang, O. V. Yazyev, S. Meng, and E. Kaxiras, Physical review letters 102, 157201 (2009).
  * Yazyev [2010] O. V. Yazyev, Reports on Progress in Physics 73, 056501 (2010).
  * Ortiz  _et al._ [2019] R. Ortiz, R. Á. Boto, N. García-Martínez, J. C. Sancho-García, M. Melle-Franco, and J. Fernández-Rossier, Nano Lett. 19, 5991 (2019).
  * Ruffieux  _et al._ [2016] P. Ruffieux, S. Wang, B. Yang, C. Sánchez-Sánchez, J. Liu, T. Dienel, L. Talirz, P. Shinde, C. A. Pignedoli, D. Passerone, _et al._ , Nature 531, 489 (2016).
  * Mishra  _et al._ [2019] S. Mishra, D. Beyer, K. Eimre, J. Liu, R. Berger, O. Groning, C. A. Pignedoli, K. Müllen, R. Fasel, X. Feng, _et al._ , Journal of the American Chemical Society 141, 10621 (2019).
  * Su  _et al._ [2019] J. Su, M. Telychko, P. Hu, G. Macam, P. Mutombo, H. Zhang, Y. Bao, F. Cheng, Z.-Q. Huang, Z. Qiu, _et al._ , Science advances 5, eaav7717 (2019).
  * Su  _et al._ [2020b] J. Su, M. Telychko, S. Song, and J. Lu, Angewandte Chemie International Edition 59, 7658 (2020b).
  * Mishra  _et al._ [2021a] S. Mishra, K. Xu, K. Eimre, H. Komber, J. Ma, C. A. Pignedoli, R. Fasel, X. Feng, and P. Ruffieux, Nanoscale 13, 1624 (2021a).
  * Mishra  _et al._ [2020] S. Mishra, D. Beyer, K. Eimre, R. Ortiz, J. Fernández-Rossier, R. Berger, O. Gröning, C. A. Pignedoli, R. Fasel, X. Feng, _et al._ , Angewandte Chemie International Edition (2020).
  * Mishra  _et al._ [2021b] S. Mishra, G. Catarina, F. Wu, R. Ortiz, D. Jacob, K. Eimre, J. Ma, C. A. Pignedoli, X. Feng, P. Ruffieux, _et al._ , Nature 598, 287 (2021b).
  * Hieulle  _et al._ [2021] J. Hieulle, S. Castro, N. Friedrich, A. Vegliante, F. R. Lara, S. Sanz, D. Rey, M. Corso, T. Frederiksen, J. I. Pascual, _et al._ , Angewandte Chemie International Edition 60, 25224 (2021).
  * Cheng  _et al._ [2022] S. Cheng, Z. Xue, C. Li, Y. Liu, L. Xiang, Y. Ke, K. Yan, S. Wang, and P. Yu, Nature Communications 13, 1705 (2022).
  * Tanaka  _et al._ [2018] K. Tanaka, Y. Yokoyama, and C. Hotta, Journal of the Physical Society of Japan 87, 023702 (2018).
  * Hu  _et al._ [2020] W.-J. Hu, S.-S. Gong, H.-H. Lai, Q. Si, and E. Dagotto, Physical Review B 101, 014421 (2020).
  * Catarina and Fernández-Rossier [2022] G. Catarina and J. Fernández-Rossier, Phys. Rev. B 105, L081116 (2022).
  * Soni  _et al._ [2022] R. Soni, N. Kaushal, C. Şen, F. A. Reboredo, A. Moreo, and E. Dagotto, New Journal of Physics 24, 073014 (2022).
  * Catarina  _et al._ [2023] G. Catarina, J. C. G. Henriques, A. Molina-Sánchez, A. T. Costa, and J. Fernández-Rossier, arXiv:2306.17153 (2023).
  * Ortiz  _et al._ [2022] R. Ortiz, G. Catarina, and J. Fernández-Rossier, 2D Materials 10, 015015 (2022).
  * Jacob  _et al._ [2021] D. Jacob, R. Ortiz, and J. Fernández-Rossier, Physical Review B 104, 075404 (2021).
  * Jacob and Fernández-Rossier [2022] D. Jacob and J. Fernández-Rossier, Physical Review B 106, 205405 (2022).
  * Creutz [1999] M. Creutz, Physical review letters 83, 2636 (1999).
  * Calzado and Malrieu [2004] C. J. Calzado and J.-P. Malrieu, Physical Review B 69, 094435 (2004).
  * Malrieu  _et al._ [2014] J. P. Malrieu, R. Caballol, C. J. Calzado, C. De Graaf, and N. Guihery, Chemical reviews 114, 429 (2014).
  * MacDonald  _et al._ [1988] A. H. MacDonald, S. Girvin, and D. t. Yoshioka, Physical Review B 37, 9753 (1988).
  * Su  _et al._ [2020c] X. Su, C. Li, Q. Du, K. Tao, S. Wang, and P. Yu, Nano letters 20, 6859 (2020c).


$\dot{p}^{ka}(t)-\lambda^{k}w^{ka}(t)=\nabla_{a}f_{1}\big{(}\bar{a}^{k}(\vartheta^{k}(t)),\bar{x}^{k}(\vartheta^{k}(t))\big{)}^{*}\big{(}\lambda^{k}(v^{kx}(t)+\theta^{kx}(t))-q^{kx}(\vartheta^{k}(t)+h_{k})\big{)},$
(8.25)
$\dot{p}^{kb}(t)-\lambda^{k}w^{kb}(t)=h_{k}\nabla_{b}f_{2}\big{(}\bar{a}^{k}(\vartheta^{k}(t)),\bar{x}^{k}(\vartheta^{k}(t))\big{)}^{*}\big{(}\lambda^{k}(v^{kx}(t)+\theta^{kx}(t))-q^{kx}(\vartheta^{k}(t)+h_{k})\big{)}$
(8.26)

hold for every $t\in(t^{k}_{j},t^{k}_{j+1})$, $j=0,\ldots,k-1$, and
$i=1,\ldots,s$.

For each $k\in I\\!\\!N$, define the vector measure $\gamma^{k}$ on $[0,T]$ by

$\displaystyle\int\limits_{B}d\,\gamma^{k}$
$\displaystyle:=\int\limits_{B}\Big{(}\sum\limits_{i=1}^{s}\eta^{k}_{i}(t)\nabla^{2}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\big{)}\big{(}\lambda^{k}(v^{kx}(t)+\theta^{kx}(t))-q^{kx}(\vartheta^{k}(t)+h_{k})\Big{)}\,dt$
(8.27)
$\displaystyle+\int\limits_{B}\Big{(}\sum\limits_{i=1}^{s}\sigma^{k}_{i}(t)\nabla
g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\Big{)}\,dt$

for any Borel subset $B\subset[0,T]$, where the vector function
$\sigma^{k}(t)=(\sigma^{k}_{1}(t),\ldots,\sigma^{k}_{s}(t))$ is taken from
(8.17). Normalizing the nontriviality condition (7.42), we get that

$\lambda^{k}+\lVert p^{k}(T)\rVert+\lVert
q^{ku}(0)\rVert+\sum\limits_{j=0}^{k}\lVert
p^{kd}_{j}\rVert=1,\,\,\,k\in\mathbb{N},$ (8.28)

which tells us that all the sequential terms in (8.28) are uniformly bounded.
Thus we have without loss of generality that there exists $\lambda\geq 0$ with
$\lambda^{k}\to\lambda$ as $k\to\infty$. It follows from (8.18) that

$\int\limits_{0}^{T}\lVert\varsigma^{k}(t)\rVert\,dt=\sum\limits_{j=0}^{k-1}\int\limits_{t^{k}_{j}}^{t^{k}_{j+1}}\lVert\varsigma^{k}(t)\rVert\,dt=\sum\limits_{j=0}^{k-1}\int\limits_{t^{k}_{j}}^{t^{k}_{j+1}}\dfrac{\lVert
p^{kd}_{j+1}\rVert}{h_{k}}\,dt=\sum\limits_{j=0}^{k-1}\lVert
p^{kd}_{j+1}\rVert=\sum\limits_{j=1}^{k}\lVert p^{kd}_{j}\rVert\leq 1,\quad
k\in I\\!\\!N.$ (8.29)

Arguing further as in the proof of [11, Theorem 6.1] gives us the following:
$\bullet$ The boundedness and the uniform bounded variations of
$\\{q^{k}(\cdot)\\}$. Hence the Helly theorem ensures the existence of a
function $q(\cdot)$ with bounded variation on $[0,T]$ such that $q^{k}(t)\to
q(t)$ as $k\to\infty$ for all $t\in[0,T]$.
$\bullet$ The boundedness of the sequence $\\{p^{k}(\cdot)\\}$ in
$W^{1,2}([0,T],\mathbb{R}^{3n+m+d})$, and hence its weak compactness in this
space. It follows therefore from the aforementioned Mazur theorem and basic
real analysis that there exists $p(\cdot)\in
W^{1,2}([0,T],\mathbb{R}^{3n+m+d})$ such that a sequence of convex
combinations of $\dot{p}^{k}(t)$ converges to $\dot{p}(t)$ for a.e.
$t\in[0,T]$. Then we arrive at (8.6) by passing to the limit in (8.22)–(8.26)
as $k\to\infty$.

Our next step is to show that the sequence $\\{\gamma^{k}(\cdot)\\}$ is
bounded in ${\cal C}([0,T],\mathbb{R}^{n})^{*}$. To proceed, take any Borel
subset $B\subset[0,T]$ and deduce from (7.32) that

$\displaystyle\Big{\|}\int\limits_{B}\Big{(}\sum\limits_{i=1}^{s}\eta^{k}_{i}(t)\nabla_{u}^{2}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\big{)}\big{(}\lambda^{k}(v^{kx}(t)+\theta^{kx}(t))-q^{kx}(\vartheta^{k}(t)+h_{k})\Big{)}\,dt$
$\displaystyle+\int\limits_{B}\Big{(}\sum\limits_{i=1}^{s}\sigma^{k}_{i}(t)\nabla_{u}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\Big{)}\,dt\Big{\|}$
$\displaystyle\leq\Big{\|}\int\limits_{0}^{T}\Big{(}\sum\limits_{i=1}^{s}\eta^{k}_{i}(t)\nabla_{u}^{2}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\big{)}\big{(}\lambda^{k}(v^{kx}(t)+\theta^{kx}(t))-q^{kx}(\vartheta^{k}(t)+h_{k})\Big{)}\,dt$
$\displaystyle+\int\limits_{0}^{T}\Big{(}\sum\limits_{i=1}^{s}\sigma^{k}_{i}(t)\nabla_{u}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\Big{)}\,dt\Big{\|}$
$\displaystyle=\Big{\|}\sum\limits_{j=0}^{k-1}\int\limits_{t_{j}}^{t_{j+1}}\Big{(}\sum\limits_{i=1}^{s}\eta^{k}_{i}(t)\nabla_{u}^{2}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\big{)}\big{(}\lambda^{k}(v^{kx}(t)+\theta^{kx}(t))-q^{kx}(\vartheta^{k}(t)+h_{k})\Big{)}\,dt$
$\displaystyle+\sum\limits_{j=0}^{k-1}\int\limits_{t_{j}}^{t_{j+1}}\Big{(}\sum\limits_{i=1}^{s}\sigma^{k}_{i}(t)\nabla_{u}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\Big{)}\,dt\Big{\|}=\Big{\|}\sum\limits_{j=0}^{k-1}\int\limits_{t_{j}}^{t_{j+1}}\Big{(}\dfrac{p^{u}_{j+1}-p^{u}_{j}}{h_{k}}-\lambda^{k}w^{u}_{j}\Big{)}\,dt\Big{\|}$
$\displaystyle=\Big{\|}\sum\limits_{j=0}^{k-1}(p^{u}_{j+1}-p^{u}_{j}-h_{k}\lambda^{k}w^{u}_{j})\Big{\|}\leq\sum\limits_{j=0}^{k-1}\|p^{u}_{j+1}-p^{u}_{j}\|+\lambda^{k}\sum\limits_{j=0}^{k-1}\|h_{k}w^{u}_{j}\|$
$\displaystyle=\sum\limits_{j=0}^{k-1}\|q^{u}(t_{j+1})-q^{u}(t_{j})\|+\lambda^{k}\sum\limits_{j=0}^{k-1}\|h_{k}w^{u}_{j}\|.$

Then from (8.28), the imposed structure of the running cost (8.1) with the
Lipschitzian functions therein, and the uniform bounded variations of
$\\{q^{u}(\cdot)\\}$ we deduce that the sequence of

$\displaystyle\int\limits_{B}\Big{(}\sum\limits_{i=1}^{s}\eta^{k}_{i}(t)\nabla_{u}^{2}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\big{)}\big{(}\lambda^{k}(v^{kx}(t)+\theta^{kx}(t))-q^{kx}(\vartheta^{k}(t)+h_{k})\Big{)}\,dt$
$\displaystyle+\int\limits_{B}\Big{(}\sum\limits_{i=1}^{s}\sigma^{k}_{i}(t)\nabla_{u}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\Big{)}\,dt$

is uniformly bounded on $[0,T]$. In the same way it follows from (7.30),
(8.28), and the uniform bounded variations of $\\{q^{x}(\cdot)\\}$ that the
sequence of

$\displaystyle\int\limits_{B}\Big{(}\sum\limits_{i=1}^{s}\eta^{k}_{i}(t)\nabla_{x}^{2}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\big{)}\big{(}\lambda^{k}(v^{kx}(t)+\theta^{kx}(t))-q^{kx}(\vartheta^{k}(t)+h_{k})\Big{)}\,dt$
$\displaystyle+\int\limits_{B}\Big{(}\sum\limits_{i=1}^{s}\sigma^{k}_{i}(t)\nabla_{x}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\Big{)}\,dt$

is also uniformly bounded on $[0,T]$. This verifies the boundedness in ${\cal
C}([0,T],\mathbb{R}^{n})^{*}$ of the sequence $\\{\gamma^{k}\\}$.

Thus we get from the weak∗ sequential compactness of the unit ball in ${\cal
C}([0,T],\mathbb{R}^{n})^{*}$ that there exists a measure $\gamma\in{\cal
C}([0,T],\mathbb{R}^{n})^{*}$ such that $\\{\gamma^{k}\\}$ weak∗ converges to
$\gamma$ along a subsequence (without relabeling). This allows us to derive
from (8.27) and the construction of the measures $\gamma^{k}$ in (8.27) that
for all $t\in[0,T]$ the following convergence holds:

$\displaystyle\int\limits_{[t,T]}\Big{(}\sum\limits_{i=1}^{s}\eta^{k}_{i}(t)\nabla^{2}_{x}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\big{)}\big{(}\lambda^{k}(v^{kx}(t)+\theta^{kx}(t))-q^{kx}(\vartheta^{k}(t)+h_{k})\Big{)}\,dt$
$\displaystyle+\int\limits_{[t,T]}\Big{(}\sum\limits_{i=1}^{s}\sigma^{k}_{i}(t)\nabla_{x}g_{i}\big{(}\bar{x}^{k}(\vartheta^{k}(t))-\bar{u}^{k}(\vartheta^{k}(t))\Big{)}\,dt\to\int\limits_{[t,T]}d\gamma(s)\;\mbox{
as }\;k\to\infty.$

Employing now (7.25), the definition of $\varsigma^{k}(\cdot)$ in (8.18), as
well as the convergence $q^{ky}(t)\to q^{y}(t)$ and $\theta^{ky}(t)\to 0$ for
a.e. $t\in[0,T]$ as $k\to\infty$ implies that

$\varsigma^{k}(t)\to q^{y}(t)\;\mbox{ a.e. }\;t\in[0,T]\;\mbox{ as
}\;k\to\infty.$

Observe by (5.3) and (7.2) that
$\dfrac{\theta^{ky}_{j}}{h_{k}}=\dfrac{1}{h_{k}}\displaystyle\int\limits_{t^{k}_{j}}^{t^{k}_{j+1}}(f_{2}(\bar{b}^{k}_{j},\bar{x}^{k}_{j})-f_{2}(\bar{b}(t),\bar{x}(t)))\,dt$,
which leads us to the estimates

$\displaystyle\lVert\varsigma^{k}(t)\rVert$ $\displaystyle\leq\lVert
q^{ky}(\vartheta^{k}(t)+h_{k})\rVert+\dfrac{\lambda^{k}}{h_{k}}\int\limits_{t^{k}_{j}}^{t^{k}_{j+1}}\lVert
f_{2}(\bar{b}^{k}_{j},\bar{x}^{k}_{j})-f_{2}(\bar{b}(t),\bar{x}(t))\rVert\,dt$
$\displaystyle\leq\lVert
q^{ky}(\vartheta^{k}(t)+h_{k})\rVert+L_{1}\dfrac{\lambda^{k}}{h_{k}}\int\limits_{t^{k}_{j}}^{t^{k}_{j+1}}\lVert\bar{x}^{k}_{j}-\bar{x}(t)\rVert\,dt+L_{2}\dfrac{\lambda^{k}}{h_{k}}\int\limits_{t^{k}_{j}}^{t^{k}_{j+1}}\lVert\bar{b}^{k}_{j}-\bar{b}(t)\rVert\,dt.$

By the optimality of $\bar{z}^{k}$ in $(P_{k})$ and the constraints in (5.5)
we get

$\lVert\varsigma^{k}(t)\rVert\leq\lVert
q^{ky}(\vartheta^{k}(t)+h_{k})\rVert+\lambda^{k}(L_{1}+L_{2})\dfrac{\varepsilon}{2}.$

The boundedness of $\\{q^{k}(\cdot)\\}$ and the normalization condition (8.28)
yield the boundedness of $\\{\vartheta^{k}(\cdot)\\}$. Then it follows from
the Lebesgue dominate convergence theorem that

$\int\limits_{[t,T]}\varsigma^{k}(s)\,ds\to\int\limits_{[t,T]}q^{y}(s)\,ds\;\mbox{
as }\;k\to\infty\;\mbox{ for all }\;t\in[0,T],$ (8.30)

and that the function $q^{y}(\cdot)$ belongs to $L^{1}([0,T],\mathbb{R}^{n})$.
Furthermore, for all $t\in[0,T]$ we have

$\displaystyle\;\quad\Big{\lVert}\int\limits_{[t,T]}\xi^{x}(\vartheta^{k}(s))\varsigma^{k}(s)\,ds-\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))q^{y}(s)\,ds\Big{\rVert}$
$\displaystyle=\Big{\lVert}\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}^{k}(\vartheta^{k}(s)),\bar{x}^{k}(\vartheta^{k}(s)))\varsigma^{k}(s)\,ds-\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))q^{y}(s)\,ds\Big{\rVert}$
$\displaystyle\leq\Big{\lVert}\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}^{k}(\vartheta^{k}(s)),\bar{x}^{k}(\vartheta^{k}(s)))\vartheta^{k}(s)\,ds-\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))\varsigma^{k}(s)\,ds\Big{\rVert}$
$\displaystyle+\Big{\lVert}\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))\varsigma^{k}(s)\,ds-\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))q^{y}(s)\,ds\Big{\rVert}.$

Note first that
$\Big{\lVert}\displaystyle\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}^{k}(\vartheta^{k}(s)),\bar{x}^{k}(\vartheta^{k}(s)))\varsigma^{k}(s)\,ds-\displaystyle\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))\varsigma^{k}(s)\,ds\Big{\rVert}\to
0$ as $k\to\infty$ due to the uniform convergence
$\bar{b}^{k}(\cdot)\to\bar{b}(\cdot)$ and
$\bar{x}^{k}(\cdot)\to\bar{x}(\cdot)$ on $[0,T]$ as $k\to\infty$, continuous
differentiability of $f_{2}(\cdot,\cdot)$, and the uniform boundedness of
$\displaystyle\int\limits_{0}^{T}\lVert\varsigma^{k}(t)\rVert\,dt$ over $k\in
I\\!\\!N$ that follows from (8.29). For the second term in the above estimate
we get
$\Big{\lVert}\displaystyle\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))\varsigma^{k}(s)\,ds-\displaystyle\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))q^{y}(s)\,ds\Big{\rVert}\to
0$ as $k\to\infty$ by the convergence of $\varsigma^{k}(\cdot)\to
q^{y}(\cdot)$ as $k\to\infty$ in $L^{1}([0,T],\mathbb{R}^{n})$. Combining this
yields

$\int\limits_{[t,T]}\xi^{x}(\vartheta^{k}(s))\varsigma^{k}(s)\,ds\to\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))q^{y}(s)\,ds\;\mbox{
as }\;k\to\infty.$

In the same way we also verify that

$\int\limits_{[t,T]}\xi^{b}(\vartheta^{k}(s))\varsigma^{k}(s)\,ds\to\int\limits_{[t,T]}\nabla_{b}f_{2}(\bar{b}(s),\bar{x}(s))q^{y}(s)\,ds\;\mbox{
as }\;k\to\infty$

and hence arrive at the adjoint relationships (8.8) by passing to the limit in
(8.20) as $k\to\infty$, where the justification of the “except a countable
subset” can be done similarly to [38, p. 325]. The passage to the limit in
(8.11) as $k\to\infty$ brings us to the integral equation

$q^{y}(t)=p^{y}(t)-\int\limits_{[t,T]}q^{y}(s)\,ds\;\mbox{ a.e. }\;t\in[0,T],$

which tells us that $q^{y}(\cdot)$ is an absolutely continuous function
satisfying

$\dot{q}^{y}(t)=\dot{p}^{y}(t)+q^{y}(t)\;\mbox{ a.e. }\;t\in[0,T].$ (8.31)

It follows from (8.6) and (8.8) that

$\dot{p}^{y}(t)=\lambda v^{x}(t)-q^{x}(t)\;\mbox{ a.e. }\;t\in[0,T],\;\mbox{
and}$ (8.32)
$q^{x}(t)=p^{x}(t)+\int\limits_{[t,T]}d\,\gamma(s)+\int\limits_{[t,T]}\nabla_{x}f_{2}(\bar{b}(s),\bar{x}(s))q^{y}(s)\,ds\;\mbox{
a.e. }\;t\in[0,T].$ (8.33)

Substituting (8.33) into (8.32) and then (8.32) into (8.31), we obtain the
extended Volterra condition (8.9). Furthermore, it follows from (7.25)–(7.28)
and the definition of $\vartheta(t)$ in (8.19) that for each $k\in I\\!\\!N$
we have

$q^{ku}(\vartheta^{k}(t)+h_{k})=\lambda^{k}(\theta^{ku}(t)+v^{ku}(t)),$ (8.34)
$q^{ka}(\vartheta^{k}(t)+h_{k})=\lambda^{k}(\theta^{ka}(t)+v^{ka}(t)),\;\mbox{
and }\;q^{kb}(\vartheta^{k}(t)+h_{k})=\lambda^{k}(\theta^{kb}(t)+v^{kb}(t)).$
(8.35)

Involving now (8.2), the assumptions on $l_{2}$, $l_{3}$, $l_{4}$, and the
Lebesgue dominated convergence theorem gives us both conditions in (8.7) by
passing to the limit in (8.34) and (8.35). The proof of the complementarity
slackness conditions (8.4), (8.5) and the endpoint conditions in (8.10)–(8.12)
is similar to the one in [11, Theorem 6.1], and we skip it for brevity. Taking
into account the convergence in (8.30) and the relationship

$\int\limits_{0}^{T}\lVert\varsigma^{k}(t)\rVert\,dt=\sum\limits_{j=0}^{k-1}\lVert
p^{kd}_{j+1}\rVert,$

we arrive at the nontriviality condition (8.13) by passing to the limit in
(8.28) as $k\to\infty$. The verification of the enhanced nontriviality
conditions in (8.14) and (8.15) under the imposed additional assumptions can
be easily proved while arguing by contraposition. This therefore completes the
proof of the theorem. $\hfill\Box$

##  9 Applications to Control of Non-Regular Electrical Circuits

This section is entirely devoted to applications of the necessary optimality
conditions for integro-differential sweeping control problems obtained in
Theorem 8.1 to controlled models that appear in non-regular electrical
circuits with ideal diodes. However, the dynamics in such models has been
described via (uncontrolled) integro-differential sweeping processes (see,
e.g., [1, 5]), we are not familiar with formulations of any optimization
and/or control problems for such systems. This is done in the first two
examples below in different frameworks. The third example of its own interest
provides a complete solution of a two-dimensional optimal control problem of
the type modeled above by using the obtained necessary optimality conditions.

Figure 1: Electrical circuit with resistors,Inductances, capacitors and ideal
diodes (RLCD). Figure 2: (RLCD) circuit with controlled voltage source.

######  Example 9.1 (optimal control of non-regular electrical circuits with
controlled current source)

Consider the electrical system depicted in Figure 1 that is composed of two
resistors $R_{1}\geq 0$ and $R_{2}\geq 0$ with voltage/current laws
$V_{R_{j}}=R_{j}x_{j}$ ($j=1,2$), two inductors $L_{1}\geq 0$ and $L_{2}\geq
0$ with voltage/current laws $V_{L_{j}}=L_{j}\dot{x}_{j}$ ($j=1,2$), two
capacitors $C_{1}>0$ and $C_{2}>0$ with voltage/current laws
$V_{C_{j}}=\frac{1}{C_{j}}\int x_{j}(t)\,dt$ $(j=1,2$), two ideal diodes with
characteristics $0\leq-V_{D_{j}}\perp i_{j}\geq 0$ and an absolutely
continuous current source $i:[0,T]\rightarrow R$. Using Kirchhoff’s laws, we
have

$\left\\{\begin{array}[]{l}V_{R_{1}}+V_{R_{2}}+V_{L_{1}}+V_{C_{1}}=-V_{D_{1}}\in-
N_{\mathbb{R}_{+}}(x_{1}-i),\\\
V_{R_{1}}-V_{R_{2}}+V_{L_{2}}+V_{C_{2}}=-V_{D_{2}}\in-
N_{\mathbb{R}_{+}}(x_{2}).\end{array}\right.$

Therefore, the dynamics of this circuit is given by

$\begin{array}[]{ll}\overbrace{\begin{pmatrix}-\dot{x}_{1}(t)\\\
-\dot{x}_{2}(t)\end{pmatrix}}^{-\dot{x}(t)}&\in
N_{C(t)}(x(t))+\overbrace{\begin{pmatrix}\frac{R_{1}+R_{2}}{L_{1}}&-\frac{R_{2}}{L_{1}}\\\
-\frac{R_{2}}{L_{2}}&\frac{R_{1}+R_{2}}{L_{2}}\end{pmatrix}}^{A_{1}}\overbrace{\begin{pmatrix}x_{1}(t)\\\
x_{2}(t)\end{pmatrix}}^{x(t)}\\\
&+\displaystyle\int\limits_{0}^{t}\bigg{[}\overbrace{\begin{pmatrix}\frac{1}{L_{1}C_{1}}&0\\\
0&\frac{1}{L_{2}C_{2}}\end{pmatrix}}^{A_{2}}\overbrace{\begin{pmatrix}x_{1}(s)\\\
x_{2}(s)\end{pmatrix}}^{x(s))}+\overbrace{\begin{pmatrix}\frac{1}{L_{1}C_{1}}i(s)\\\
0\end{pmatrix}}^{b(s)}\bigg{]}\,ds.\end{array}$ (9.1)

Put $u(t):=(i(t),0)^{*}$, $b(t)=(\frac{1}{L_{1}C_{1}}i(t),0)^{*}$ indicating
vector columns and denote $C(t)=C(u(t)):=u(t)+[0,\infty)\times[0,\infty)$,
$f_{1}(x):=A_{1}x$, $f_{2}(b,x):=A_{2}x+b$, and $x(0):=(i(0),0)^{*}$. Thus
(9.1) can be rewritten in the form of the controlled sweeping dynamics (1.2)
as

$-\dot{x}(t)\in
N_{C(u(t))}(x(t))+f_{1}(x(t))+\displaystyle\int\limits_{0}^{t}f_{2}(b(s),x(s))\,ds\;\mbox{
a.e. }\;t\in[0,T],\;x(0)=x_{0}\in C(0)$

with control functions acting in the moving set and the integral part of the
dynamics. The cost functional in this problem is naturally formulated by:
minimize

$J[x,u,b]:=\frac{\lambda_{T}}{2}(x_{1}(T)-i(0))^{2}+\frac{\lambda_{P}}{2}\displaystyle\int\limits_{0}^{T}(x_{1}(t)-i(0))^{2}\,dt+\frac{\lambda_{I}}{2}\displaystyle\int\limits_{0}^{T}b_{1}^{2}(t)\,dt,$

where $\lambda_{T}$, $\lambda_{Q}$, and $\lambda_{I}$ are nonnegative
constants not equal to zero simultaneously. The necessary optimality
conditions obtained in Theorem 8.1 can be readily applied to this class of
optimal control problems.

The next model describes a class of dynamic processes for non-regular electric
circuits with an input controlled voltage source, which is dual of the current
source considered in Example 9.1. In this model we have control actions
entering the dynamic perturbation term of the sweeping process.

######  Example 9.2 (optimal control of non-regular electrical circuit with
controlled voltage source)

Consider the electrical system shown in Figure 2 that is composed of a
resistor $R\geq 0$ with voltage/current law $V_{R}=R(x_{3}-x_{2})$, two
inductors $L_{1}\geq 0$, $L_{2}\geq 0$ with voltage/current laws
$V_{L_{j}}=L_{j}\dot{x}_{j}$ ($j=1,2$), a capacitor $C>0$ with voltage/current
law $V_{C}=\dfrac{1}{C}\displaystyle\int(x_{2}(t)-x_{1}(t))\,dt$, two ideal
diodes with characteristics $0\leq-V_{D_{j}}\perp x_{j}\geq 0$, and two
voltage sources $v_{j}(t)$ ($j=1,2$). Employing Kirchhoff’s laws tells us that

$\left\\{\begin{array}[]{l}V_{L_{1}}-V_{C}-v_{1}(t)=-V_{D_{1}}\in-
N_{\mathbb{R}_{+}}(x_{1}),\\\ V_{L_{2}}+V_{C}-V_{R}=-V_{D_{2}}\in-
N_{\mathbb{R}_{+}}(x_{2}),\\\ V_{R}-v_{2}(t)=0.\end{array}\right.$

In the system above, we substitute the third equation into the second one and
get

$-\begin{pmatrix}\dot{x}_{1}(t)\\\ \dot{x}_{2}(t)\end{pmatrix}\in
N_{\mathbb{R}^{2}_{+}}\begin{pmatrix}x_{1}(t)\\\
x_{2}(t)\end{pmatrix}+\overbrace{\begin{pmatrix}-\frac{1}{L_{1}}&0\\\
0&-\frac{1}{L_{2}}\end{pmatrix}}^{A_{1}}\begin{pmatrix}v_{1}(t)\\\
v_{2}(t)\end{pmatrix}+\int\limits_{0}^{t}\overbrace{\begin{pmatrix}\frac{1}{L_{1}C}&-\frac{1}{L_{1}C}\\\
-\frac{1}{L_{2}C}&\frac{1}{L_{2}C}\end{pmatrix}}^{A_{2}}\begin{pmatrix}x_{1}(s)\\\
x_{2}(s)\end{pmatrix}\,ds.$

The goal is to start from a given state $(x_{1}(0),x_{2}(0))$ and reach the
other state $(x_{1}(T),x_{2}(T))$ with minimizing the input energy. Denoting
$C:=\\{x=(x_{1},x_{2})\in\mathbb{R}^{2}\;|\;g_{1}(x)=x_{1}\geq
0,\,\,\,g_{2}(x)=x_{2}\geq 0\\}$, $u(t)=(0,0)$, $x_{0}\in C$,
$a(t)=(a_{1}(t),a_{2}(t))$ with $(a_{1}(t),a_{2}(t)):=(v_{1}(t),v_{2}(t)$,
$f_{1}(a,x):=A_{1}a$, and $f_{2}(b,x):=A_{2}x$, the formulated problem can be
written in the form of our basic problem $(P)$ as follows:

$\mbox{minimize
}\;J[x,a]:=\dfrac{x_{2}^{2}(T)}{2}+\frac{1}{2}\int\limits_{0}^{T}[a_{1}^{2}(t)+a_{2}^{2}(t)]\,dt$
(9.2)

subject to the controlled integro-differential sweeping process

$-\dot{x}(t)\in
N_{C}(x(t))+f_{1}(a(t),x(t))+\displaystyle\int\limits_{0}^{t}f_{2}(x(s))\,ds\;\mbox{
a.e. }\;[0,T],\;x(0)=x_{0}\in C.$ (9.3)

It follows from the existence result of Theorem 3.2 that the optimal control
problem formulated in (9.2) and (9.3) admits an optimal solution. Furthermore,
the structure of the cost functional in (9.2) together with Proposition 3.1
yields the uniqueness of the optimal pair $(\bar{x}(\cdot),\bar{a}(\cdot))$ in
problem (9.2), (9.3). This allows us to find the optimal solution to the
formulated optimal control problems by using the necessary optimality
conditions Theorem 8.1, provided that a feasible solution determined by these
conditions is unique.

Finally, we present a numerical example that demonstrates how to use the
necessary optimality conditions of Theorem 8.1 to solve a particular case of
the optimal control problem formulated in Example 9.2 with the given data of
controlled integro-differential sweeping inclusion. The corresponding example
for the first model described in this section can be constructed similarly,
while we skip it for brevity.

######  Example 9.3 (calculating optimal solutions of two-dimensional sweeping
control model)

Consider the optimal control problem defined in (9.2) and (9.3) with the
following data:

$\begin{array}[]{ll}T=1,\;x_{0}=(1,1),\;g(x)=(g_{1}(x),g_{2}(x)):=(x_{1},x_{2}),\;l_{0}(x,a):=\dfrac{a_{1}^{2}+a_{2}^{2}}{2}.\;\varphi(x)=\dfrac{x^{2}_{2}}{2},\\\
A_{1}:=\begin{pmatrix}-1&0\\\ 0&-1\end{pmatrix},\;\mbox{ and
}\;A_{2}=\begin{pmatrix}1&-1\\\ -1&1\end{pmatrix}.\end{array}$ (9.4)

We seek for solutions to problem (9.4) such that

$\bar{x}(t)>0\;\mbox{ for all }\;t\in[0,1)\;\mbox{ and }\;\bar{x}(1)\in{\rm
bd}\,C.$ (9.5)

Depending on whether the diodes are blocking (off) or conducting (on),
condition (9.5) provides the following three possibilities that are listed
below as modes:
$\bullet$ Mode 1: For all $t\in[0,1)$ both diodes are on, i.e., $x_{1}>0$
$(V_{D_{1}}=0)$ and $x_{2}>0$ $(V_{D_{2}}=0)$. Furthermore, at the ending time
$T=1$ the first diode is off while the second one is on, i.e., $x_{1}(1)=0$
and $x_{2}(1)>0$ $(V_{D_{2}}=0)$.
$\bullet$ Mode 2: For all $t\in[0,1)$ both diodes are on, i.e., $x_{1}>0$
$(V_{D_{1}}=0)$ and $x_{2}>0$ $(V_{D_{2}}=0)$, while at the ending time $T=1$
the first diode is on and the second one is off, i.e., $x_{1}(1)>0$
$(V_{D_{1}}=0)$ and $x_{2}(1)=0.$
$\bullet$ Mode 3: For all $t\in[0,1)$ both diodes are on, i.e., $x_{1}>0$
$(V_{D_{1}}=0)$ and $x_{2}>0$ $(V_{D_{2}}=0)$, while at the ending time $T=1$
both diodes are off, i.e., $x_{1}(1)=0$ and $x_{2}(1)=0.$

Applying the necessary optimality conditions of Theorem 8.1 gives us a number
$\lambda\geq 0$, functions $\eta_{i}(\cdot)\in L^{2}([0,1],\mathbb{R}_{+})$ as
$i=1,2$ well-defined at $t=1$, triples
$p(\cdot)=(p^{x}(\cdot),p^{y}(\cdot),p^{a}(\cdot))\in
W^{1,2}([0,T];\mathbb{R}^{6})$ and
$q(\cdot)=(q^{x}(\cdot),q^{y}(\cdot),q^{a}(\cdot))$ with values in
$\mathbb{R}^{6}$ and of bounded variations on $[0,T]$, as well as a pair
$w(\cdot)=(w^{x}(\cdot),w^{a}(\cdot))\in L^{2}([0,T],\mathbb{R}^{4})$ such
that

  1. 1.

$w(t)=(0,0,\bar{a}(t)),\,\;v(t)=(0,0,0)\quad\text{a.e.}\quad t\in[0,1].$

  2. 2.

$\begin{cases}\dot{\bar{x}}_{1}(t)-\bar{a}_{1}(t)+\displaystyle\int\limits_{0}^{t}(\bar{x}_{1}(s)-\bar{x}_{2}(s))\,ds=\eta_{1}(t)\;\mbox{
a.e. }\;t\in[0,1],&\\\
\dot{\bar{x}}_{2}(t)-\bar{a}_{2}(t)-\displaystyle\int\limits_{0}^{t}(\bar{x}_{1}(s)-\bar{x}_{2}(s))\,ds=\eta_{2}(t)\;\mbox{
a.e. }\;t\in[0,1].\end{cases}$

  3. 3.

$\bar{x}_{i}(t)>0\Longrightarrow\eta_{i}(t)=0\;\mbox{ a.e.
}\;t\in[0,1],\;\;i=1,2.$

  4. 4.

$\eta_{i}(t)>0\Longrightarrow q_{i}^{x}(t)=0\;\mbox{ a.e.
}\;t\in[0,1),\;i=1,2$.

  5. 5.

$(\dot{p}^{x}(t),\dot{p}^{y}(t),\dot{p}^{a}(t))=((0,0),-q^{x}(t),\lambda\bar{a}(t)-A_{1}q^{x}(t))\;\mbox{
a.e. }\;t\in[0,1]$.

  6. 6.

$q^{a}(t)=(0,0)\;\mbox{ a.e. }\;t\in[0,1]$.

  7. 7.

$(q^{x}(t),q^{y}(t),q^{a}(t))=(p^{x}(t),p^{y}(t),p^{a}(t))-\displaystyle\int\limits_{[t,1]}(-d\gamma(s)-A_{2}q^{y}(s)\,ds,q^{y}(s)\,ds,(0,0))\;\mbox{
a.e. }\;t\in[0,1]$.

  8. 8.

$\begin{cases}\dot{q}_{1}^{y}(t)=-p_{1}^{x}(t)-\displaystyle\int\limits_{[t,1]}d\gamma_{1}(s)-\displaystyle\int\limits_{[t,1]}(q_{1}^{y}(s)-q_{2}^{y}(s))\,ds+q_{1}^{y}(t)\;\mbox{
a.e. }\;t\in[0,1],\\\
\dot{q}_{2}^{y}(t)=-p_{2}^{x}(t)-\displaystyle\int\limits_{[t,1]}d\gamma_{2}(s)+\displaystyle\int\limits_{[t,1]}(q_{1}^{y}(s)-q_{2}^{y}(s))\,ds+q_{2}^{y}(t)\;\mbox{
a.e. }\;t\in[0,1].\end{cases}$

  9. 9.

$-p_{1}^{x}(1)+\eta_{1}(1)=0,\,\,-p_{2}^{x}(1)+\eta_{2}(1)=\lambda\bar{x}_{2}(1)$.

  10. 10.

$-\eta(1)\in N_{C}(\bar{x}(1)).$

  11. 11.

$p^{y}(1)=p^{a}(1)=(0,0)$.

  12. 12.

$\lambda+|\eta_{1}(1)|+|\eta_{2}(1)|+\displaystyle\int\limits^{1}_{0}\|q^{y}(t)\|\,dt>0$.

It follows directly from items 5–7 that

$p^{x}(t)\equiv p^{x}(1)\;\mbox{ and }\;\lambda\bar{a}(t)=-q^{x}(t)\;\mbox{
for all }\;t\in[0,1].$

Proceeding now similarly to [10, Example 1] gives us
$\displaystyle\int\limits_{[t,1]}d\,\gamma(s)=\gamma(\\{1\\})$. Then

$q^{x}(t)=p^{x}(1)+\gamma(\\{1\\})+\displaystyle\int\limits_{[t,1]}A_{2}q^{y}(s)\,ds,\quad
t\in[0,1].$ (9.6)

Solving the integro-differential system in item 8 with $q^{y}(1)=(0,0)$, we
get on $[0,1]$ that

$\begin{cases}\displaystyle\int\limits_{[t,1]}q_{1}^{y}(s)\,ds=\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})+p^{x}_{2}(1)+\gamma_{2}(\\{1\\})}{2}e^{t-1}+\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{6}e^{-(t-1)}\\\
+\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{12}e^{2t-2}-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{4}-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})+p^{x}_{2}(1)+\gamma_{2}(\\{1\\})}{2}t,\\\
\\\
\displaystyle\int\limits_{[t,1]}q_{2}^{y}(s)\,ds=\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})+p^{x}_{2}(1)+\gamma_{2}(\\{1\\})}{2}e^{t-1}-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{6}e^{-(t-1)}\\\
-\displaystyle\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{12}e^{2t-2}+\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{4}-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})+p^{x}_{2}(1)+\gamma_{2}(\\{1\\})}{2}t.\end{cases}$

Then it follows from the above relationships, equation (9.6), and the
definition of $A_{2}$ that

$\begin{cases}q_{1}^{x}(t)=p^{x}_{1}(1)+\gamma_{1}(\\{1\\})+\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{3}e^{-(t-1)}\vspace*{0.3cm}\\\
+\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{6}e^{2t-2}-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{2}&\vspace*{0.5cm}\textstyle\\\
q_{2}^{x}(t)=p^{x}_{2}(1)+\gamma_{2}(\\{1\\})-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{3}e^{-(t-1)}\vspace*{0.3cm}\\\
-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{6}e^{2t-2}+\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{2}.\end{cases}$

This leads us to the representations

$\begin{cases}\lambda\bar{a}_{1}(t)=-q_{1}^{x}(t)=-p^{x}_{1}(1)-\gamma_{1}(\\{1\\})-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{3}e^{-(t-1)}\vspace*{0.3cm}\\\
-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{6}e^{2t-2}+\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{2},&\vspace*{0.5cm}\textstyle\\\
\lambda\bar{a}_{2}(t)=-q_{2}^{x}(t)=-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})+\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{3}e^{-(t-1)}\vspace*{0.3cm}\\\
+\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{6}e^{2t-2}-\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})-p^{x}_{2}(1)-\gamma_{2}(\\{1\\})}{2},\end{cases}$

for all $t\in[0,1]$. Suppose now for simplicity that $\lambda>0$ and denote
$v_{1}:=\dfrac{p^{x}_{1}(1)+\gamma_{1}(\\{1\\})}{\lambda}$,
$v_{2}:=\dfrac{p^{x}_{2}(1)+\gamma_{2}(\\{1\\})}{\lambda}$. Then the
nontriviality condition in item 12 holds automatically, and we have

$\begin{cases}\bar{a}_{1}(t)=-v_{1}-\dfrac{v_{1}-v_{2}}{3}e^{-(t-1)}-\dfrac{v_{1}-v_{2}}{6}e^{2t-2}+\dfrac{v_{1}-v_{2}}{2},&\vspace*{0.5cm}\textstyle\\\
\bar{a}_{2}(t)=-v_{2}+\dfrac{v_{1}-v_{2}}{3}e^{-(t-1)}+\dfrac{v_{1}-v_{2}}{6}e^{2t-2}-\dfrac{v_{1}-v_{2}}{2}.\end{cases}$
(9.7)

Substituting these expressions into item 2 and then using item 3 yield

$\begin{cases}\dot{\bar{x}}_{1}(t)+v_{1}+\dfrac{v_{1}-v_{2}}{3}e^{-(t-1)}+\dfrac{v_{1}-v_{2}}{6}e^{2t-2}-\dfrac{v_{1}-v_{2}}{2}+\displaystyle\int\limits_{0}^{t}(\bar{x}_{1}(s)-\bar{x}_{2}(s))\,ds=0,&\vspace*{0.5cm}\textstyle\\\
\dot{\bar{x}}_{2}(t)+v_{2}-\dfrac{v_{1}-v_{2}}{3}e^{-(t-1)}-\dfrac{v_{1}-v_{2}}{6}e^{2t-2}+\dfrac{v_{1}-v_{2}}{2}-\displaystyle\int\limits_{0}^{t}(\bar{x}_{1}(s)-\bar{x}_{2}(s))\,ds=0.\end{cases}$

Solving the obtained integro-differential system with $\bar{x}(0)=(1,1)$ gives
us

$\begin{cases}\bar{x}_{1}(t)=1-\dfrac{v_{1}+v_{2}}{2}t-\dfrac{v_{1}-v_{2}}{18}e^{2t-2}+\dfrac{v_{1}-v_{2}}{9}e^{1-t}+\dfrac{1}{3}\cos(\sqrt{2}t)\Big{(}\dfrac{v_{1}-v_{2}}{6}e^{-2}-\dfrac{v_{1}-v_{2}}{3}e\Big{)}\vspace*{0.3cm}\\\
-\dfrac{1}{\sqrt{2}}\sin(\sqrt{2}t)\Big{(}\dfrac{v_{1}-v_{2}}{18}e^{-2}+\dfrac{2(v_{1}-v_{2})}{9}e\Big{)},&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{2}(t)=1-\dfrac{v_{1}+v_{2}}{2}t+\dfrac{v_{1}-v_{2}}{18}e^{2t-2}-\dfrac{v_{1}-v_{2}}{9}e^{1-t}-\dfrac{1}{3}\cos(\sqrt{2}t)\Big{(}\dfrac{v_{1}-v_{2}}{6}e^{-2}-\dfrac{v_{1}-v_{2}}{3}e\Big{)}\vspace*{0.3cm}\\\
+\dfrac{1}{\sqrt{2}}\sin(\sqrt{2}t)\Big{(}\dfrac{v_{1}-v_{2}}{18}e^{-2}+\dfrac{2(v_{1}-v_{2})}{9}e\Big{)},&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{3}(t)=\bar{a}_{2}(t)+\bar{x}_{2}(t).\end{cases}$ (9.8)

By the second condition in (9.5) we have the following three possibilities :

  1. (i)

If $\bar{x}_{1}(1)=0$, then $v_{1}=v_{2}-\dfrac{v_{2}-1}{c}$, where
$c:=\dfrac{4}{9}-\dfrac{\cos(\sqrt{2})}{9}\Big{(}\dfrac{e^{-2}}{2}-e\Big{)}+\dfrac{\sin(\sqrt{2})}{9\sqrt{2}}\Big{(}\dfrac{e^{-2}}{2}+2e\Big{)}$.
It is easy to deduce from (9.7) and (9.8) that

$\begin{cases}\bar{a}_{1}(t)=-v_{2}+\dfrac{v_{2}-1}{3c}e^{-(t-1)}+\dfrac{v_{2}-1}{6c}e^{2t-2}+\dfrac{v_{2}-1}{2c},&\vspace*{0.5cm}\textstyle\\\
\bar{a}_{2}(t)=-v_{2}-\dfrac{v_{2}-1}{3c}e^{-(t-1)}-\dfrac{v_{2}-1}{6c}e^{2t-2}+\dfrac{v_{2}-1}{2c}.\end{cases}$
$\begin{cases}\bar{x}_{1}(t)=1-v_{2}t+\dfrac{v_{2}-1}{c}\Big{[}\dfrac{t}{2}+\dfrac{e^{2t-2}}{18}-\dfrac{e^{1-t}}{9}-\dfrac{1}{9}\cos(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}-e\Big{)}+\dfrac{1}{9\sqrt{2}}\sin(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}+2e\Big{)}\Big{]},&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{2}(t)=1-v_{2}t+\dfrac{v_{2}-1}{c}\Big{[}\dfrac{t}{2}-\dfrac{e^{2t-2}}{18}+\dfrac{e^{1-t}}{9}+\dfrac{1}{9}\cos(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}-e\Big{)}-\dfrac{1}{9\sqrt{2}}\sin(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}+2e\Big{)}\Big{]},&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{3}(t)=\bar{a}_{2}(t)+\bar{x}_{2}(t).\end{cases}$

Observe further that
$\bar{x}_{2}(t)+\bar{x}_{1}(t)=2-2v_{2}t+\dfrac{v_{2}-1}{c}t$. Then we have
that $\bar{x}_{2}(1)=2-2v_{2}+\dfrac{v_{2}-1}{c}$, and that the cost
functional in (9.2) reduces to

$\displaystyle J[\bar{x},\bar{a}]$
$\displaystyle=\dfrac{1}{2}\Big{(}2-2v_{2}+\dfrac{v_{2}-1}{c}\Big{)}^{2}+\dfrac{1}{2}\int\limits_{0}^{1}\Big{[}-v_{2}+\dfrac{v_{2}-1}{c}\Big{(}\dfrac{1}{2}+\dfrac{e^{-(t-1)}}{3}+\dfrac{e^{2t-2}}{6}\Big{)}\Big{]}^{2}\,dt$
$\displaystyle+\dfrac{1}{2}\int\limits_{0}^{1}\Big{[}-v_{2}+\dfrac{v_{2}-1}{c}\Big{(}\dfrac{1}{2}-\dfrac{e^{-(t-1)}}{3}-\dfrac{e^{2t-2}}{6}\Big{)}\Big{]}^{2}\,dt.$

This function clearly achieves its absolute minimum at $v_{2}=0.5988481275$.
Hence

$\begin{cases}\bar{a}_{1}(t)=-0.5988481275-\dfrac{0.4011518725}{3c}e^{-(t-1)}-\dfrac{0.4011518725}{6c}e^{2t-2}-\dfrac{0.4011518725}{2c},&\vspace*{0.5cm}\textstyle\\\
\bar{a}_{2}(t)=-0.5988481275+\dfrac{0.4011518725}{3c}e^{-(t-1)}+\dfrac{0.4011518725}{6c}^{2t-2}-\dfrac{0.4011518725}{2c},\end{cases}$
$\begin{cases}\bar{x}_{1}(t)=1-0.5988481275.t-\dfrac{0.4011518725}{c}\Big{[}\dfrac{t}{2}+\dfrac{e^{2t-2}}{18}-\dfrac{e^{1-t}}{9}-\dfrac{1}{9}\cos(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}-e\Big{)}\vspace*{0.3cm}\\\
+\dfrac{1}{9\sqrt{2}}\sin(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}+2e\Big{)}\Big{]},&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{2}(t)=1-0.5988481275.t-\dfrac{0.4011518725}{c}\Big{[}\dfrac{t}{2}-\dfrac{e^{2t-2}}{18}+\dfrac{e^{1-t}}{9}+\dfrac{1}{9}\cos(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}-e\Big{)}\vspace*{0.3cm}\\\
-\dfrac{1}{9\sqrt{2}}\sin(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}+2e\Big{)}\Big{]},&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{3}(t)=\bar{a}_{2}(t)+\bar{x}_{2}(t)\;\mbox{ for all
}\;t\in[0,1].\end{cases}$

  2. (ii)

If $\bar{x}_{2}(1)=0$, then $v_{1}=v_{2}+\dfrac{v_{2}-1}{c}$, where
$c:=-\dfrac{5}{9}-\dfrac{\cos(\sqrt{2})}{9}\Big{(}\dfrac{e^{-2}}{2}-e\Big{)}+\dfrac{\sin(\sqrt{2})}{9\sqrt{2}}\Big{(}\dfrac{e^{-2}}{2}+2e\Big{)}$.
Arguing similar to the above Case (i), we get that

$\begin{cases}\bar{a}_{1}(t)=-1.056787399-\dfrac{0.056787399}{3c}e^{-(t-1)}-\dfrac{0.056787399}{6c}e^{2t-2}-\dfrac{0.056787399}{2c},&\vspace*{0.5cm}\textstyle\\\
\bar{a}_{2}(t)=-1.056787399+\dfrac{0.056787399}{3c}e^{-(t-1)}+\dfrac{0.056787399}{6c}e^{2t-2}-\dfrac{0.056787399}{2c},\end{cases}$
$\begin{cases}\bar{x}_{1}(t)=1-1.056787399.t+\dfrac{0.056787399}{c}\Big{[}-\dfrac{t}{2}-\dfrac{e^{2t-2}}{18}+\dfrac{e^{1-t}}{9}+\dfrac{1}{9}\cos(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}-e\Big{)}\vspace*{0.3cm}\\\
-\dfrac{1}{9\sqrt{2}}\sin(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}+2e\Big{)}\Big{]},&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{2}(t)=1-1.056787399.t+\dfrac{0.056787399}{c}\Big{[}-\dfrac{t}{2}+\dfrac{e^{2t-2}}{18}-\dfrac{e^{1-t}}{9}-\dfrac{1}{9}\cos(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}-e\Big{)}\vspace*{0.3cm}\\\
+\dfrac{1}{9\sqrt{2}}\sin(\sqrt{2}t)\Big{(}\dfrac{e^{-2}}{2}+2e\Big{)}\Big{]},&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{3}(t)=\bar{a}_{2}(t)+\bar{x}_{2}(t)\;\mbox{ for all
}\;t\in[0,1].\end{cases}$

  3. (iii)

Let $\bar{x}_{1}(1)=0$ and $\bar{x}_{2}(1)=0$. Then it follows from Case (i)
that $v_{1}=v_{2}-\dfrac{v_{2}-1}{c}$ and $0=2-2v_{2}+\dfrac{v_{2}-1}{c}$,
where
$c:=\dfrac{4}{9}-\dfrac{\cos(\sqrt{2})}{9}\Big{(}\dfrac{e^{-2}}{2}-e\Big{)}+\dfrac{\sin(\sqrt{2})}{9\sqrt{2}}\Big{(}\dfrac{e^{-2}}{2}+2e\Big{)}$.
Hence $v_{1}=v_{2}=1$.
Furthermore, we deduce from (9.7) and (9.8) the following expressions valid
for all $t\in[0,1]$:

$\begin{cases}\bar{a}_{1}(t)=\bar{a}_{2}(t)=-1,&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{1}(t)=\bar{x}_{2}(t)=1-t,&\vspace*{0.5cm}\textstyle\\\
\bar{x}_{3}(t)=\bar{a}_{2}(t)+\bar{x}_{2}(t)=-t.\end{cases}$

Thus the obtained necessary optimality conditions from Theorem 8.1 with
$\lambda\neq 0$ allows us to fully compute the unique optimal solution of the
integro-differential sweeping control problem formulated in (9.2) and (9.3)
with the initial data given in (9.4) such that condition (9.5) is satisfied.

## References

  * [1] V. Acary, O. Bonnefon and B. Brogliato, Nonsmooth Modeling and Simulation for Switched Circuits, Springer, Dordrecht, The Netherlands, 2011.
  * [2] L. Adam and J.V. Outrata, On optimal control of a sweeping process coupled with an ordinary differential equation, Discrete Contin. Dyn. Syst. Ser. B 19 (2014), 2709–2738.
  * [3] S. Adly, F. Nacry and L. Thibault, Preservation of prox-regularity of sets with applications to constrained optimization, SIAM J. Optim. 26 (2016), 448–473.
  * [4] C.E. Arroud and G. Colombo, A maximum principle of the controlled sweeping process, Set-Valued Var. Anal. 26 (2018), 607–629.
  * [5] A. Bouach, T. Haddad and L. Thibault, Nonconvex integro-differential sweeping process with applications, arXiv:2102.11987 (2021).
  * [6] Y. Brenier. W. Gangbo, G. Savare and M. Westdickenberg, Sticky particle dynamics with interactions, J. Math. Pures Appl. 99 (2013), 577–617.
  * [7] B. Brogliato and A. Tanwani, Dynamical systems coupled with monotone set-valued operators: formalisms, applications, well-posedness, and stability, SIAM Rev. 62 (2020), 3–129.
  * [8] M. Brokate and P. Krejčí, Optimal control of ODE systems involving a rate independent variational inequality, Discrete Contin. Dyn. Syst. Ser. B 18 (2013), 331–348.
  * [9] T.H. Cao, G. Colombo, B.S. Mordukhovich and D. Nguyen, Optimization of fully controlled sweeping processes, to appear in J. Diff. Eqs.; arXiv:2012.11722v1.
  * [10] T.H. Cao and B.S. Mordukhovich, Optimality conditions for a controlled sweeping process with applications to the crowd motion model, Discrete Contin. Dyn. Syst. Ser. B 21 (2017), 267–306.
  * [11] T.H. Cao and B.S. Mordukhovich, Optimal control of a nonconvex perturbed sweeping process, J. Diff. Eqs. 266 (2019), 1003–1050.
  * [12] T.H. Cao and B.S. Mordukhovich, Applications of optimal control of a nonconvex sweeping process to optimization of the planar crowd motion model, Discrete Contin. Dyn. Syst. Ser. B 24 (2019), 4191–4216.
  * [13] D.C. Clark, A short proof of discrete Gronwall inequality, Disc. Appl. Math. 16 (1987), 279–281.
  * [14] F.H. Clarke, Yu.S. Ledyaev, R.J. Stern and P.R. Wolenski, Nonsmooth Analysis and Control Theory, Springer, New York, 1998.
  * [15] G. Colombo, R. Henrion, N.D. Hoang and B.S. Mordukhovich, Optimal control of the sweeping process, Dyn. Contin. Discrete Impuls. Syst. Ser. B 19 (2012), 117–159.
  * [16] G. Colombo, R. Henrion, N.D. Hoang and B.S. Mordukhovich, Optimal control of the sweeping process over polyhedral controlled sets, J. Diff. Eqs. 260 (2016), 3397–3447.
  * [17] G. Colombo and C. Kozaily, Existence and uniqueness of solutions for an integral perturbation of Moreau’s sweeping process, J. Convex Anal. 27 (2020), 227–236.
  * [18] G. Colombo, B.S. Mordukhovich and D. Nguyen, Optimization of a perturbed sweeping process by discontinuous controls, SIAM J. Control Optim. 58 (2020), 2678–2709.
  * [19] G. Colombo and L. Thibault, Prox-regular sets and applications, in: D.Y. Gao and D. Motreanu (Eds.), Handbook of Nonconvex Analysis, International Press, Boston, 2010, pp. 99–182.
  * [20] M.d.R. de Pinho, M.M.A. Ferreira and G.V. Smirnov, Optimal control involving sweeping processes, Set-Valued Var. Anal. 27 (2019), 523–548.
  * [21] T. Donchev, E. Farkhi and B.S. Mordukhovich, Discrete approximations, relaxation, and optimization of one-sided Lipschitzian differential inclusions in Hilbert spaces, J. Diff. Eqs. 243 (2007), 301–328.
  * [22] J. F. Edmond and L. Thibault, Relaxation of an optimal control problem involving a perturbed sweeping process, Math. Program. 104 (2005), 347–373.
  * [23] R. Henrion, J.V. Outrata and T. Surowiec, On the coderivative of normal cone mappings to inequality systems, Nonlinear Analysis 71 (2009), 1213–1226.
  * [24] N.D. Hoang and B.S. Mordukhovich, Extended Euler-Lagrange and Hamiltonian formalisms in optimal control of sweeping processes with controlled sweeping sets, J. Optim. Theory Appl. 180 (2019), 256–289.
  * [25] M. Kunze and M.D.P. Monteiro Marques, An introduction to Moreau‘s sweeping process, in: B. Brogliato (Ed.), Impacts in Mechanical Systems, in: Lecture Notes in Phys., vol. 551, 2000, Springer, Berlin, 2000, pp. 1–60.
  * [26] B.S. Mordukhovich, Sensitivity analysis in nonsmooth optimization, in: D.A. Field and V. Komkov (Eds.), Theoretical Aspects of Industrial Design, in SIAM Proc. Applied Math., vol. 58, 1992, pp. 32–46.
  * [27] B.S. Mordukhovich, Discrete approximations and refined Euler-Lagrange conditions for differential inclusions, SIAM J. Control Optim. 33 (1995), 882–915.
  * [28] B.S. Mordukhovich, Variational Analysis and Generalized Differentiation, I: Basic Theory, Springer, Berlin, 2006.
  * [29] B.S. Mordukhovich, Variational Analysis and Generalized Differentiation, II: Applications, Springer, Berlin, 2006.
  * [30] B.S. Mordukhovich, Variational Analysis and Applications, Springer, Cham, Switzerland, 2018.
  * [31] B.S. Mordukhovich, Optimal control of Lipschitzian and discontinuous differential inclusions with various applications, Proc. Inst. Math. Mech. Azer. Acad. Sci. 45 (2019), 52–74.
  * [32] B. S. Mordukhovich and D. Nguyen, Discrete approximations and optimal control of nonsmooth perturbed sweeping processes, J. Convex Anal. 28, No. 2, 2021.
  * [33] J.J. Moreau, On unilateral constraints, friction and plasticity, in: G. Capriz and G. Stampacchia (Eds.), New Variational Techniques in Mathematical Physics, Proceedings of C.I.M.E. Summer Schools, Cremonese, Rome, 1974, pp. 173–322.
  * [34] R.T. Rockafellar, Maximal monotone relations and the second derivatives in nonsmootrh analysis, Ann. Inst. H. Poincaré: Analyse Non Linéare 2, 167–184.
  * [35] R.T. Rockafellar and R.J-B. Wets, Variational Analysis, Springer, Berlin, 1998.
  * [36] A.A. Tolstonogov, Control sweeping process, J. Convex Anal. 23 (2016), 1099–1123.
  * [37] J. Venel, A numerical scheme for a class of sweeping processes, Numer. Math. 118 (2011), 451–484.
  * [38] R.B. Vinter, Optimal Control, Birkhaüser, Boston, 2000.
  * [39] V. Zeidan, C. Nour and H. Saoud, A nonsmooth maximum principle for a controlled nonconvex sweeping process, J. Diff. Eqs. 269 (2020), 9531–9582.