# A short note on the appearance of the simplest
antilinear ODE in several physical contexts

Dmitry Ponomarev1,2,3

###### Abstract

In this short note, we review several one-dimensional problems such as those
involving linear Schrödinger equation, variable-coefficient Helmholtz
equation, Zakharov-Shabat system and Kubelka-Munk equations. We show that they
all can be reduced to solving one simple antilinear ordinary differential
equation $u^{\prime}\left(x\right)=f\left(x\right)\overline{u\left(x\right)}$
or its nonhomogeneous version
$u^{\prime}\left(x\right)=f\left(x\right)\overline{u\left(x\right)}+g\left(x\right)$,
$x\in\left(0,x_{0}\right)\subset\mathbb{R}$. We point out some of the
advantages of the proposed reformulation and call for further investigation of
the obtained ODE.

##  1 Introduction

11footnotetext: Institute of Analysis and Scientific Computing, Vienna
University of Technology (TU Wien), Austria22footnotetext: St. Petersburg
Department of Steklov Mathematical Institute of Russian Academy of Sciences,
Russia33footnotetext: Contact<EMAIL_ADDRESS>
Many physical phenomena can be directly described by or reduced to systems of
differential equations having certain structural properties. Restricting
ourselves here to linear one-dimensional settings, we are concerned with a
pair of first-order ODEs whose matrix is antidiagonal with complex-conjugate
elements. Namely, given $x_{0}\in\mathbb{R}$ and complex-valued function
$f\left(x\right)$, we consider the equation

$U^{\prime}\left(x\right)=\left(\begin{array}[]{cc}0&f\left(x\right)\\\
\overline{f\left(x\right)}&0\end{array}\right)U\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right),$ (1)

as well as its nonhomogeneous analog

$U^{\prime}\left(x\right)=\left(\begin{array}[]{cc}0&f\left(x\right)\\\
\overline{f\left(x\right)}&0\end{array}\right)U\left(x\right)+G\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right)\text{,}$ (2)

where
$U\left(x\right)\equiv\left(u_{1}\left(x\right),u_{2}\left(x\right)\right)^{T}\in\mathbb{C}^{2}$
is an unknown solution-vector,
$G\left(x\right)\equiv\left(g_{1}\left(x\right),g_{2}\left(x\right)\right)^{T}\in\mathbb{C}^{2}$
is a given vector-function, and each of equations (1)–(2) is supplemented by
the initial condition $U\left(0\right)=U_{0}\in\mathbb{C}^{2}$. Here and
onwards, we employ the notation $\overline{\cdot}$ to denote complex
conjugation.

Similarly to Hamiltonian, Dirac and more general canonical systems (see e.g.
[12]), equations (1)–(2) constitute an important class of dynamical systems
for two reasons. On the one hand, as we shall further see, formulations of
several important problems are reducible to either (1) or (2). On the other
hand, these systems are close to being exactly solvable in the following
sense. Let us focus on (1) and consider the more general system

$U^{\prime}\left(x\right)=\left(\begin{array}[]{cc}p\left(x\right)&r\left(x\right)\\\
s\left(x\right)&q\left(x\right)\end{array}\right)U\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right).$ (3)

We note that the diagonal elements in the matrix of (3) can be removed by the
exponential multiplier transform. Namely, by setting

$V\left(x\right):=\left(\begin{array}[]{cc}e^{-\int_{0}^{x}p\left(\tau\right)\textrm{d}\tau}&0\\\
0&e^{-\int_{0}^{x}q\left(\tau\right)\textrm{d}\tau}\end{array}\right)U\left(x\right),$

one can observe that $V\left(x\right)$ satisfies

$V^{\prime}\left(x\right)=\left(\begin{array}[]{cc}0&r\left(x\right)\exp\left(-\int_{0}^{x}\left[p\left(\tau\right)-q\left(\tau\right)\right]\textrm{d}\tau\right)\\\
s\left(x\right)\exp\left(\int_{0}^{x}\left[p\left(\tau\right)-q\left(\tau\right)\right]\textrm{d}\tau\right)&0\end{array}\right)V\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right),$ (4)

with the initial condition $V\left(0\right)=U_{0}$. Now, if the anti-diagonal
elements of the matrix in the right-hand side of (4) are equal, i.e.

$r\left(x\right)\exp\left(-\int_{0}^{x}\left[p\left(\tau\right)-q\left(\tau\right)\right]\textrm{d}\tau\right)=s\left(x\right)\exp\left(\int_{0}^{x}\left[p\left(\tau\right)-q\left(\tau\right)\right]\textrm{d}\tau\right)=:c_{1}\left(x\right),$
(5)

then the solution can be written explicitly as

$V\left(x\right)=\left[\cosh\left(\int_{0}^{x}c_{1}\left(\tau\right)\textrm{d}\tau\right)I+\sinh\left(\int_{0}^{x}c_{1}\left(\tau\right)\textrm{d}\tau\right)S\right]U_{0},$

where

$I:=\left(\begin{array}[]{cc}1&0\\\ 0&1\end{array}\right),\hskip
10.00002pt\hskip 10.00002ptS:=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right).$ (6)

However, condition (5), which amounts to the assumption

$s\left(x\right)=r\left(x\right)\exp\left(-2\int_{0}^{x}\left[p\left(\tau\right)-q\left(\tau\right)\right]\textrm{d}\tau\right),$
(7)

may be too restrictive. Indeed, in view of multiple possible similarity
transformations allowing to rewrite (3) in different equivalent forms, we want
to have a clearly identifiable matrix structure which should be, on the one
hand, immediately recognisable and, on the other hand, leading to a solution
simplification or even an explicit solution. Such an identifiable structure
may be, for example, a pairwise relation between some of the elements of the
matrix in (3). The explicit solvability condition, nevertheless, plays against
any visible structural property of the matrix: even though assumption (7)
leads to a closed-form solution, it implies a rather complicated relation
between the matrix elements. Condition (7) is very specific and thus unlikely
to be satisfied for any easily describable class of matrix elements unless, of
course, $p\equiv q$, which would then also entail that $r\equiv s$. This
exactly solvable case with equal diagonal and anti-diagonal elements in (3)
may be sometimes valuable but it does not seem to be the one that covers many
fruitful applications.

It turns out that condition (7) has an analog which is less stringent in form
of the matrix elements in (3), with more pertinence to important physical
contexts, and, at the same time, it still leads to a significant
simplification of the solution procedure (and, at least in some cases, also to
closed-form solutions). This condition reads

$s\left(x\right)=\overline{r\left(x\right)}\exp\left(-2\text{
Re}\int_{0}^{x}\left[p\left(\tau\right)-q\left(\tau\right)\right]\textrm{d}\tau\right).$
(8)

Despite the similarity to (7), condition (8) is easier to satisfy while
preserving a visible matrix structure. Indeed, if
$p\left(x\right)-q\left(x\right)$ is a purely imaginary function (e.g., in
particular, when $p\equiv\overline{q}$), the complicated exponential factor in
(8) disappears. In this case, the implied condition
$s\left(x\right)=\overline{r\left(x\right)}$ is clearly indentifiable but far
from being trivial since, as we shall see, it covers a variety of different
practical applications. This reasoning motivates us to consider (1) as well as
its nonhogeneous analog (2).

The plan of this note is as follows. Section 2 is dedicated to the
transformation of (1) and (2) into the homogeneous antilinear ODE

$u^{\prime}\left(x\right)=f\left(x\right)\overline{u\left(x\right)},\hskip
10.00002ptx\in\left(0,x_{0}\right),$ (9)

and its nonhomogeneous analog

$u^{\prime}\left(x\right)=f\left(x\right)\overline{u\left(x\right)}+g\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right),$ (10)

respectively. Next, in Section 3, we outline some relevant physical
applications, i.e. problems which, upon appropriate transformations, can be
recast as (1) or (2) and are thus reducible to formulations involving
antilinear ODE (9), or, in one case, its nonhomogeneous version (10). Finally,
in Section 4, we conclude with some remarks on how antilinear ODEs can be
constructively addressed further and briefly mention a couple of other
applications.

##  2 Transformation of an antidiagonal problem into an antilinear ODE

###  2.1 Homogeneous case: from (1) to (9)

We consider (1) supplemented with the initial data
$U\left(0\right)=U_{0}\equiv\left(u_{1}^{0},u_{2}^{0}\right)^{T}$, where
$u_{1}^{0}$, $u_{2}^{0}\in\mathbb{C}$ are arbitrary constants, and we devise a
transformation that allows construction of the solution of system (1) in terms
of the solution of an antilinear ODE of the form (9).

Let us first motivate our approach to construction of such a transformation.
To this effect, given a complex-valued function $f\left(x\right)$, it is
instructive to consider the elementary differential equation for
$v\left(x\right)$

$v^{\prime}\left(x\right)=f\left(x\right)v\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right),$ (11)

with the initial condition $v\left(0\right)=1$. On the one hand, (11) is in
separable form and hence can be integrated directly to yield the solution

$v\left(x\right)=\exp\left(\int_{0}^{x}f\left(\tau\right)\textrm{d}\tau\right).$
(12)

On the other hand, rewriting (11) in the integral form

$v\left(x\right)=1+\int_{0}^{x}f\left(\tau\right)v\left(\tau\right)\textrm{d}\tau,\hskip
10.00002ptx\in\left(0,x_{0}\right),$

the Picard iterative process gives

$v\left(x\right)=1+\int_{0}^{x}f\left(\tau\right)\textrm{d}\tau+\int_{0}^{x}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}f\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}+\int_{0}^{x}f\left(\tau_{3}\right)\int_{0}^{\tau_{3}}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}f\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}\textrm{d}\tau_{3}+\ldots.$
(13)

Comparison of (12) with (13) results in important identities

$\exp\left(\int_{0}^{x}f\left(\tau\right)\textrm{d}\tau\right)=1+\int_{0}^{x}f\left(\tau\right)\textrm{d}\tau+\int_{0}^{x}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}f\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}+\int_{0}^{x}f\left(\tau_{3}\right)\int_{0}^{\tau_{3}}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}f\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}\textrm{d}\tau_{3}+\ldots,$
$\sinh\left(\int_{0}^{x}f\left(\tau\right)\textrm{d}\tau\right)=\int_{0}^{x}f\left(\tau\right)\textrm{d}\tau+\int_{0}^{x}f\left(\tau_{3}\right)\int_{0}^{\tau_{3}}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}f\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}\textrm{d}\tau_{3}+\ldots,$
(14)
$\cosh\left(\int_{0}^{x}f\left(\tau\right)\textrm{d}\tau\right)=1+\int_{0}^{x}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}f\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}+\int_{0}^{x}f\left(\tau_{4}\right)\int_{0}^{\tau_{4}}f\left(\tau_{3}\right)\int_{0}^{\tau_{3}}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}f\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}\textrm{d}\tau_{3}\textrm{d}\tau_{4}+\ldots,$
(15)

where we used the identity
$\exp\left(z\right)=\cosh\left(z\right)+\sinh\left(z\right)$,
$z\in\mathbb{C}$, and the parity argument to split the terms: $\sinh$ is an
odd function and hence (14) may contain only odd number of multiplicative
instances of $f$, and similarly, (15) may contain only terms with even number
of multiplications by $f$ due to $\cosh$ being an even function.

Now, similarly to (14)–(15), let us consider the following quantities

$S_{f}\left(x\right):=\int_{0}^{x}f\left(\tau_{1}\right)\textrm{d}\tau_{1}+\int_{0}^{x}f\left(\tau_{3}\right)\int_{0}^{\tau_{3}}\overline{f\left(\tau_{2}\right)}\int_{0}^{\tau_{2}}f\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}\textrm{d}\tau_{3}+\ldots,$
(16)
$C_{f}\left(x\right):=1+\int_{0}^{x}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}\overline{f\left(\tau_{1}\right)}\textrm{d}\tau_{1}\textrm{d}\tau_{2}+\int_{0}^{x}f\left(\tau_{4}\right)\int_{0}^{\tau_{4}}\overline{f\left(\tau_{3}\right)}\int_{0}^{\tau_{3}}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}\overline{f\left(\tau_{1}\right)}\textrm{d}\tau_{1}\textrm{d}\tau_{2}\textrm{d}\tau_{3}\textrm{d}\tau_{4}+\ldots.$
(17)

Let us show that (16)–(17) are inherent to an algebraic structure underlying
(1). To this effect, we rewrite (1) in the integral form

$U\left(x\right)=U_{0}+\int_{0}^{x}A\left(\tau\right)U\left(\tau\right)\textrm{d}\tau,\hskip
10.00002ptx\in\left(0,x_{0}\right),\hskip 10.00002pt\hskip
10.00002ptA\left(\tau\right):=\left(\begin{array}[]{cc}0&f\left(\tau\right)\\\
\overline{f\left(\tau\right)}&0\end{array}\right),$ (18)

and note that

$A\left(\tau_{2}\right)A\left(\tau_{1}\right)=\left(\begin{array}[]{cc}f\left(\tau_{2}\right)\overline{f\left(\tau_{1}\right)}&0\\\
0&\overline{f\left(\tau_{2}\right)}f\left(\tau_{1}\right)\end{array}\right),\hskip
10.00002ptA\left(\tau_{3}\right)A\left(\tau_{2}\right)A\left(\tau_{1}\right)=\left(\begin{array}[]{cc}0&f\left(\tau_{3}\right)\overline{f\left(\tau_{2}\right)}f\left(\tau_{1}\right)\\\
\overline{f\left(\tau_{3}\right)}f\left(\tau_{2}\right)\overline{f\left(\tau_{1}\right)}&0\end{array}\right),$
$A\left(\tau_{4}\right)A\left(\tau_{3}\right)A\left(\tau_{2}\right)A\left(\tau_{1}\right)=\left(\begin{array}[]{cc}f\left(\tau_{4}\right)\overline{f\left(\tau_{3}\right)}f\left(\tau_{2}\right)\overline{f\left(\tau_{1}\right)}&0\\\
0&\overline{f\left(\tau_{4}\right)}f\left(\tau_{3}\right)\overline{f\left(\tau_{2}\right)}f\left(\tau_{1}\right)\end{array}\right),\hskip
10.00002pt\ldots.$

Therefore, writing out Picard iterations for solving (18), we obtain

$U\left(x\right)=\left(\begin{array}[]{cc}C_{f}\left(x\right)&S_{f}\left(x\right)\\\
S_{\overline{f}}\left(x\right)&C_{\overline{f}}\left(x\right)\end{array}\right)U_{0}=\left(\begin{array}[]{cc}C_{f}\left(x\right)&S_{f}\left(x\right)\\\
\overline{S_{f}\left(x\right)}&\overline{C_{f}\left(x\right)}\end{array}\right)U_{0}.$
(19)

Furthermore, it is easy to see from (16)–(17) that $S_{f}\left(x\right)$,
$C_{f}\left(x\right)$ obey the following intertwining relation

$C_{f}^{\prime}\left(x\right)=f\left(x\right)\overline{S_{f}\left(x\right)},\hskip
10.00002pt\hskip
10.00002ptS_{f}^{\prime}\left(x\right)=f\left(x\right)\overline{C_{f}\left(x\right)},\hskip
10.00002pt\hskip 10.00002ptx\in\left(0,x_{0}\right),$ (20)

and the conditions $S_{f}\left(0\right)=0$, $C_{f}\left(0\right)=1$.
Introducing another pair of functions

$Z_{+}\left(x\right):=C_{f}\left(x\right)+S_{f}\left(x\right),\hskip
10.00002pt\hskip
10.00002ptZ_{-}\left(x\right):=C_{f}\left(x\right)-S_{f}\left(x\right),$

we decouple (20) as

$Z_{+}^{\prime}\left(x\right)=f\left(x\right)\overline{Z_{+}\left(x\right)},\hskip
10.00002ptx\in\left(0,x_{0}\right),\hskip 10.00002pt\hskip
10.00002ptZ_{+}\left(0\right)=1,$ (21)
$Z_{-}^{\prime}\left(x\right)=-f\left(x\right)\overline{Z_{-}\left(x\right)},\hskip
10.00002ptx\in\left(0,x_{0}\right),\hskip 10.00002pt\hskip
10.00002ptZ_{-}\left(0\right)=1.$ (22)

Equations (21)–(22) are two separate instances of the initial-value problem
featuring the antilinear ODE given by (9). Solution of this ODE would thus
yield the solutions of (21)–(22) and, consequently, also of (20), providing
$S_{f}\left(x\right)$, $C_{f}\left(x\right)$ appearing in (19) which furnishes
the solution of (1).

###  2.2 Nonhomogeneous case: from (2) to (10)

Let us now consider (2) with
$G\left(x\right)\equiv\left(g_{1}\left(x\right),g_{2}\left(x\right)\right)^{T}$
and subject to the initial condition
$U\left(0\right)=U_{0}\equiv\left(u_{1}^{0},u_{2}^{0}\right)^{T}$. We are
going to show that, in particular case where

$g_{2}\left(x\right)=i\,\overline{g_{1}\left(x\right)},\hskip 10.00002pt\hskip
10.00002ptu_{2}^{0}=i\,\overline{u_{1}^{0}},$ (23)

the solution of (2) can be constucted in terms of solutions of two instances
of problem (10). As we shall see in Subsection 3.2, assumption (23) will be
satisfied in at least one important practical context.

Similarly to (16)–(17), let us introduce

$\mathrm{S}_{f,h}\left(x\right):=\int_{0}^{x}f\left(\tau_{1}\right)h\left(\tau_{1}\right)\textrm{d}\tau_{1}+\int_{0}^{x}f\left(\tau_{3}\right)\int_{0}^{\tau_{3}}\overline{f\left(\tau_{2}\right)}\int_{0}^{\tau_{2}}f\left(\tau_{1}\right)h\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}\textrm{d}\tau_{3}+\ldots,$
(24) $\displaystyle\mathrm{C}_{f,h}\left(x\right):=$ $\displaystyle
h\left(x\right)+\int_{0}^{x}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}\overline{f\left(\tau_{1}\right)}h\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}$
(25)
$\displaystyle+\int_{0}^{x}f\left(\tau_{4}\right)\int_{0}^{\tau_{4}}\overline{f\left(\tau_{3}\right)}\int_{0}^{\tau_{3}}f\left(\tau_{2}\right)\int_{0}^{\tau_{2}}\overline{f\left(\tau_{1}\right)}h\left(\tau_{1}\right)\textrm{d}\tau_{1}\textrm{d}\tau_{2}\textrm{d}\tau_{3}\textrm{d}\tau_{4}+\ldots.$

Rewriting (2) in the integral form

$U\left(x\right)=U_{0}+\int_{0}^{x}G\left(\tau\right)\textrm{d}\tau+\int_{0}^{x}A\left(\tau\right)U\left(\tau\right)\textrm{d}\tau,\hskip
10.00002ptx\in\left(0,x_{0}\right),\hskip 10.00002pt\hskip
10.00002ptA\left(\tau\right):=\left(\begin{array}[]{cc}0&f\left(\tau\right)\\\
\overline{f\left(\tau\right)}&0\end{array}\right),$ (26)

it is straightforward to see that Picard iterations give

$U\left(x\right)=\left(\begin{array}[]{c}\mathrm{S}_{f,h_{2}}\left(x\right)+\mathrm{C}_{f,h_{1}}\left(x\right)\\\
\mathrm{S}_{\overline{f},h_{1}}\left(x\right)+\mathrm{C}_{\overline{f},h_{2}}\left(x\right)\end{array}\right),$

where $\mathrm{S}_{f,h}\left(x\right)$, $\mathrm{C}_{f,h}\left(x\right)$ are
as defined by (24)–(25), and

$h_{1}\left(x\right):=u_{1}^{0}+\int_{0}^{x}g_{1}\left(\tau\right)\textrm{d}\tau,\hskip
10.00002pt\hskip
10.00002pth_{2}\left(x\right):=u_{2}^{0}+\int_{0}^{x}g_{2}\left(\tau\right)\textrm{d}\tau.$
(27)

By means of differentiation of $\mathrm{S}_{f,h_{2}}\left(x\right)$ and
$\mathrm{C}_{f,h_{1}}\left(x\right)$, we obtain the following intertwining
relation

$\mathrm{C}_{f,h_{1}}^{\prime}\left(x\right)=h_{1}^{\prime}\left(x\right)+f\left(x\right)\mathrm{S}_{\overline{f},h_{1}}\left(x\right),\hskip
10.00002pt\mathrm{S}_{f,h_{2}}^{\prime}\left(x\right)=f\left(x\right)\mathrm{C}_{\overline{f},h_{2}}\left(x\right),\hskip
10.00002pt\hskip 10.00002ptx\in\left(0,x_{0}\right),$ (28)

which is to be supplemented by the conditions
$\mathrm{C}_{f,h_{1}}\left(0\right)=h_{1}\left(0\right)$,
$\mathrm{S}_{f,h_{2}}\left(0\right)=0$.

Note that, from (24)–(25),
$\mathrm{S}_{\overline{f},h_{1}}\left(x\right)=\overline{\mathrm{S}_{f,\overline{h_{1}}}\left(x\right)}$
and
$\mathrm{C}_{\overline{f},h_{2}}\left(x\right)=\overline{\mathrm{C}_{f,\overline{h_{2}}}\left(x\right)}$.
Moreover, assumption (23) entails that
$h_{2}\left(x\right)=i\,\overline{h_{1}\left(x\right)}$, and by linearity in
$h$ of $\mathrm{S}_{f,h}\left(x\right)$, $\mathrm{C}_{f,h}\left(x\right)$, we
have
$\mathrm{S}_{f,h_{2}}\left(x\right)=i\,\mathrm{S}_{f,\overline{h_{1}}}\left(x\right)$
and
$\overline{\mathrm{C}_{f,\overline{h_{2}}}\left(x\right)}=i\,\overline{\mathrm{C}_{f,h_{1}}\left(x\right)}$.
Consequently, relation (28) becomes

$\mathrm{C}_{f,h_{1}}^{\prime}\left(x\right)=h_{1}^{\prime}\left(x\right)+f\left(x\right)\overline{\mathrm{S}_{f,\overline{h_{1}}}\left(x\right)},\hskip
10.00002pt\mathrm{S}_{f,\overline{h_{1}}}^{\prime}\left(x\right)=f\left(x\right)\overline{\mathrm{C}_{f,h_{1}}\left(x\right)},\hskip
10.00002pt\hskip 10.00002ptx\in\left(0,x_{0}\right).$ (29)

Setting

$Z_{+}\left(x\right):=\mathrm{C}_{f,h_{1}}\left(x\right)+\mathrm{S}_{f,\overline{h_{1}}}\left(x\right),\hskip
10.00002ptZ_{-}\left(x\right):=\mathrm{C}_{f,h_{1}}\left(x\right)-\mathrm{S}_{f,\overline{h_{1}}}\left(x\right),$
(30)

we obtain from (29) two decoupled ODE problems

$Z_{+}^{\prime}\left(x\right)=f\left(x\right)\overline{Z_{+}\left(x\right)}+g_{1}\left(x\right),\hskip
10.00002pt\hskip 10.00002ptZ_{+}\left(0\right)=u_{1}^{0},$ (31)
$Z_{-}^{\prime}\left(x\right)=-f\left(x\right)\overline{Z_{-}\left(x\right)}+g_{1}\left(x\right),\hskip
10.00002pt\hskip 10.00002ptZ_{-}\left(0\right)=u_{1}^{0},$ (32)

each of them are of the form (10).

##  3 Some physical contexts leading to (1) and (2)

###  3.1 Linear Schrödinger equation

Consider the stationary linear Schrödinger equation in $1D$, with a potential
$a\left(x\right)>0$,

$u^{\prime\prime}\left(x\right)+a\left(x\right)u\left(x\right)=0,\hskip
10.00002ptx\in\left(0,x_{0}\right).$ (33)

We focus here on the initial-value problem, i.e. we supplement (33) with the
boundary conditions $u\left(0\right)=u_{0}$, $u^{\prime}\left(0\right)=u_{1}$,
but boundary-value problems on $\left(0,x_{0}\right)$, with $x_{0}$ being
finite or infinite, could also be treated. We assume $a\in
C^{1}\left(\left[0,x_{0}\right]\right)$.

Introducing the vector-function

$U\left(x\right):=\left(\begin{array}[]{c}u\left(x\right)\\\
\frac{1}{a^{1/2}\left(x\right)}u^{\prime}\left(x\right)\end{array}\right),$
(34)

we observe that $U\left(x\right)$ satisfies

$U^{\prime}\left(x\right)=\left(\begin{array}[]{c}u^{\prime}\left(x\right)\\\
-\frac{a^{\prime}\left(x\right)}{2a^{3/2}\left(x\right)}u^{\prime}\left(x\right)+\frac{1}{a^{1/2}\left(x\right)}u^{\prime\prime}\left(x\right)\end{array}\right)=A\left(x\right)U\left(x\right),$
(35)

with

$A\left(x\right):=\left(\begin{array}[]{cc}0&a^{1/2}\left(x\right)\\\
-a^{1/2}\left(x\right)&-\frac{a^{\prime}\left(x\right)}{2a\left(x\right)}\end{array}\right),$

and $U\left(0\right)=\left(u_{0},u_{1}/a^{1/2}\left(0\right)\right)^{T}$.
Here, in the second inequality of (35), we used (33) to eliminate
$u^{\prime\prime}\left(x\right)$.

By writing,

$A\left(x\right)=a^{1/2}\left(x\right)\left(\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right)+\frac{a^{\prime}\left(x\right)}{2a\left(x\right)}\left(\begin{array}[]{cc}0&0\\\
0&-1\end{array}\right),$

we note that the first matrix in the right-hand side is diagonalisable as
follows

$P\left(\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right)P^{-1}=\left(\begin{array}[]{cc}i&0\\\
0&-i\end{array}\right),\hskip 10.00002pt\hskip
10.00002ptP:=\frac{1}{2^{1/2}}\left(\begin{array}[]{cc}i&1\\\
1&i\end{array}\right),\hskip
10.00002ptP^{-1}=\frac{1}{2^{1/2}}\left(\begin{array}[]{cc}-i&1\\\
1&-i\end{array}\right).$ (36)

Consequently, introducing $V\left(x\right):=PU\left(x\right)$, we multiply the
both sides of (35) by $P$ and thus transform it into

$V^{\prime}\left(x\right)=PA\left(x\right)P^{-1}V\left(x\right)=B\left(x\right)V\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right),$ (37)

with

$B\left(x\right):=a^{1/2}\left(x\right)\left(\begin{array}[]{cc}i&0\\\
0&-i\end{array}\right)-\frac{a^{\prime}\left(x\right)}{4a\left(x\right)}\left(\begin{array}[]{cc}1&-i\\\
i&1\end{array}\right)=\left(\begin{array}[]{cc}b_{0}\left(x\right)&0\\\
0&\overline{b_{0}\left(x\right)}\end{array}\right)-\frac{a^{\prime}\left(x\right)}{4a\left(x\right)}\left(\begin{array}[]{cc}0&-i\\\
i&0\end{array}\right),$
$b_{0}\left(x\right):=ia^{1/2}\left(x\right)-\frac{a^{\prime}\left(x\right)}{4a\left(x\right)},$
(38)

and supplemented by the initial condition

$V\left(0\right)=PU\left(0\right)=\frac{1}{2^{1/2}}\left(\begin{array}[]{c}iu_{0}+u_{1}/a^{1/2}\left(0\right)\\\
u_{0}+iu_{1}/a^{1/2}\left(0\right)\end{array}\right).$

Furthermore, introducing

$W\left(x\right):=\left(\begin{array}[]{cc}\exp\left(-\int_{0}^{x}b_{0}\left(\tau\right)\textrm{d}\tau\right)&0\\\
0&\exp\left(-\int_{0}^{x}\overline{b_{0}\left(\tau\right)}\textrm{d}\tau\right)\end{array}\right)V\left(x\right),$

we have

$\frac{\textrm{d}}{\textrm{d}x}\left(\begin{array}[]{cc}\exp\left(-\int_{0}^{x}b_{0}\left(\tau\right)\textrm{d}\tau\right)&0\\\
0&\exp\left(-\int_{0}^{x}\overline{b_{0}\left(\tau\right)}\textrm{d}\tau\right)\end{array}\right)=\left(\begin{array}[]{cc}-b_{0}\left(x\right)\exp\left(-\int_{0}^{x}b_{0}\left(\tau\right)\textrm{d}\tau\right)&0\\\
0&-\overline{b_{0}\left(x\right)}\exp\left(-\int_{0}^{x}\overline{b_{0}\left(\tau\right)}\textrm{d}\tau\right)\end{array}\right).$

Therefore, (37) entails

$W^{\prime}\left(x\right)=-\frac{a^{\prime}\left(x\right)}{4a\left(x\right)}\left(\begin{array}[]{cc}0&-i\exp\left(-\int_{0}^{x}\left[b_{0}\left(\tau\right)-\overline{b_{0}\left(\tau\right)}\right]\textrm{d}\tau\right)\\\
i\exp\left(\int_{0}^{x}\left[b_{0}\left(\tau\right)-\overline{b_{0}\left(\tau\right)}\right]\textrm{d}\tau\right)&0\end{array}\right)W\left(x\right),$

which, recalling (38), we can rewrite as

$W^{\prime}\left(x\right)=C\left(x\right)W\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right),$ (39)

with

$C\left(x\right):=\left(\begin{array}[]{cc}0&c_{0}\left(x\right)\\\
\overline{c_{0}\left(x\right)}&0\end{array}\right),\hskip 10.00002pt\hskip
10.00002ptc_{0}\left(x\right):=\frac{ia^{\prime}\left(x\right)}{4a\left(x\right)}\exp\left(-2i\int_{0}^{x}a^{1/2}\left(\tau\right)\textrm{d}\tau\right),$

and the initial condition

$W\left(0\right)=V\left(0\right)=\frac{1}{2^{1/2}}\left(\begin{array}[]{c}iu_{0}+u_{1}/a^{1/2}\left(0\right)\\\
u_{0}+iu_{1}/a^{1/2}\left(0\right)\end{array}\right).$

The steps described above draw from [10] (see also [4]) and provide one way to
rewrite the linear Schrödinger equation in the form (1), but this approach is
not the only one. Alternative reduction procedures may be more cumbersome but
more beneficial in practice, depending on a final goal. For instance, in [2],
the initial vectorisation of (33) is different from (34) yet other steps of
the transformation are ideologically similar.

###  3.2 Helmholtz equation

Stationary problems for the wave propagation in heterogeneous media are
described by the Helmholtz equation whose 1$D$ version is given by

$\left(\alpha\left(x\right)u^{\prime}\left(x\right)\right)^{\prime}+\beta\left(x\right)u\left(x\right)=f\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right).$ (40)

Here, $\alpha\left(x\right)$, $\beta\left(x\right)>0$ are material parameters
and $f\left(x\right)$ is the source term. As in Subsection 3.1, we suppose
that (40) is supplemented by the initial conditions $u\left(0\right)=u_{0}$,
$u^{\prime}\left(0\right)=u_{1}$. Furthermore, we assume that
$f\left(x\right)$, $u_{0}$, $u_{1}$ are all real-valued. This assumption does
not reduce generality since (40) is linear with real-valued
$\alpha$$\left(x\right)$, $\beta\left(x\right)$ and hence a real-valued
problem can be solved separately for real and imaginary parts of the solution
of the original equation.

Setting

$U\left(x\right):=\left(\begin{array}[]{c}\alpha^{1/2}\left(x\right)\beta^{1/2}\left(x\right)u\left(x\right)\\\
\alpha\left(x\right)u^{\prime}\left(x\right)\end{array}\right),$

we recast (40) in the vector form

$U^{\prime}\left(x\right)=A\left(x\right)U\left(x\right)+F_{0}\left(x\right),\hskip
10.00002ptx\in\left(0,x_{0}\right),$ (41)

with

$A\left(x\right):=\left(\begin{array}[]{cc}\frac{\left[\alpha\left(x\right)\beta\left(x\right)\right]^{\prime}}{2\alpha\left(x\right)\beta\left(x\right)}&\left(\frac{\beta\left(x\right)}{\alpha\left(x\right)}\right)^{1/2}\\\
-\left(\frac{\beta\left(x\right)}{\alpha\left(x\right)}\right)^{1/2}&0\end{array}\right),\hskip
10.00002pt\hskip 10.00002ptF_{0}\left(x\right):=\left(\begin{array}[]{c}0\\\
f\left(x\right)\end{array}\right),$

and the initial condition

$U\left(0\right)=\left(\begin{array}[]{c}\alpha^{1/2}\left(0\right)\beta^{1/2}\left(0\right)u_{0}\\\
\alpha\left(0\right)u_{1}\end{array}\right).$

We now follow the reduction steps similar to those in Subsection 3.1. We write

$A\left(x\right)=\left(\frac{\beta\left(x\right)}{\alpha\left(x\right)}\right)^{1/2}\left(\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right)+\frac{\left[\alpha\left(x\right)\beta\left(x\right)\right]^{\prime}}{2\alpha\left(x\right)\beta\left(x\right)}\left(\begin{array}[]{cc}1&0\\\
0&0\end{array}\right),$

and note that we can diagonalise the first matrix with the help of the
auxiliary constant matrix $P$ introduced in (36). Denoting
$V\left(x\right):=PU\left(x\right)$, we hence have, from (41),

$V^{\prime}\left(x\right)=B\left(x\right)V\left(x\right)+F_{1}\left(x\right),$
(42)

with

$B\left(x\right):=\left(\begin{array}[]{cc}b_{1}\left(x\right)&0\\\
0&\overline{b_{1}\left(x\right)}\end{array}\right)+\frac{\left[\alpha\left(x\right)\beta\left(x\right)\right]^{\prime}}{4\alpha\left(x\right)\beta\left(x\right)}\left(\begin{array}[]{cc}0&i\\\
-i&0\end{array}\right),\hskip
10.00002ptF_{1}\left(x\right):=PF\left(x\right)=\left(\begin{array}[]{c}\frac{f\left(x\right)}{2^{1/2}}\\\
i\frac{f\left(x\right)}{2^{1/2}}\end{array}\right),$
$b_{1}\left(x\right):=i\left(\frac{\beta\left(x\right)}{\alpha\left(x\right)}\right)^{1/2}+\frac{\left[\alpha\left(x\right)\beta\left(x\right)\right]^{\prime}}{4\alpha\left(x\right)\beta\left(x\right)},$

and the initial condition

$V\left(0\right)=PU\left(0\right)=\frac{1}{2^{1/2}}\left(\begin{array}[]{c}i\alpha^{1/2}\left(0\right)\beta^{1/2}\left(0\right)u_{0}+\alpha\left(0\right)u_{1}\\\
\alpha^{1/2}\left(0\right)\beta^{1/2}\left(0\right)u_{0}+i\alpha\left(0\right)u_{1}\end{array}\right).$

Introducing

$W\left(x\right):=\left(\begin{array}[]{cc}\exp\left(-\int_{0}^{x}b_{1}\left(\tau\right)\textrm{d}\tau\right)&0\\\
0&\exp\left(-\int_{0}^{x}\overline{b_{1}\left(\tau\right)}\textrm{d}\tau\right)\end{array}\right)V\left(x\right),$

equation (42) transforms into

$W^{\prime}\left(x\right)=C\left(x\right)W\left(x\right)+G\left(x\right),$
(43)

where

$C\left(x\right):=\left(\begin{array}[]{cc}0&c_{1}\left(x\right)\\\
\overline{c_{1}\left(x\right)}&0\end{array}\right),\hskip 10.00002pt\hskip
10.00002ptc_{1}\left(x\right):=\frac{i\left[\alpha\left(x\right)\beta\left(x\right)\right]^{\prime}}{4\alpha\left(x\right)\beta\left(x\right)}\exp\left(-2i\int_{0}^{x}\left(\frac{\beta\left(\tau\right)}{\alpha\left(\tau\right)}\right)^{1/2}\textrm{d}\tau\right),$
$G\left(x\right):=\left(\begin{array}[]{c}g_{1}\left(x\right)\\\
i\overline{g_{1}\left(x\right)}\end{array}\right),\hskip 10.00002pt\hskip
10.00002ptg_{1}\left(x\right):=\frac{f\left(x\right)}{2^{1/2}}\exp\left(-\int_{0}^{x}\left[i\left(\frac{\beta\left(\tau\right)}{\alpha\left(\tau\right)}\right)^{1/2}+\frac{\left[\alpha\left(x\right)\beta\left(x\right)\right]^{\prime}}{4\alpha\left(x\right)\beta\left(x\right)}\right]\textrm{d}\tau\right),$

and the initial condition

$W\left(0\right)=\left(\begin{array}[]{c}w_{1}^{0}\\\
i\overline{w_{1}^{0}}\end{array}\right),\hskip
10.00002ptw_{1}^{0}:=\frac{1}{2^{1/2}}\left(i\alpha^{1/2}\left(0\right)\beta^{1/2}\left(0\right)u_{0}+\alpha\left(0\right)u_{1}\right).$

Here, in relating the first and the second components of the vector
$G\left(x\right)$, and similarly $W\left(0\right)$, we employed the real-
valuedness of $\alpha\left(x\right)$, $\beta\left(x\right)$, $u_{0}$, $u_{1}$
that was discussed in the beginning of this Subsection.

It remains to observe that system (43) is such that the matrix
$C\left(x\right)$ and the vectors $G\left(x\right)$, $W\left(0\right)$ fit the
assumptions discussed in Subsection 2.2.

###  3.3 Zakharov-Shabat system

It is well-known that solution of a spectral problem with the linear
Schrödinger equation appears as an intermediate step in solving the
Korteweg–de Vries (KdV) equation by using the inverse scattering transform.
Zakharov-Shabat systems play the same role in the integrability of other
nonlinear equations [1, p.10]. In particular, the Zakharov-Shabat system

$\partial_{x}\left(\begin{array}[]{c}v_{1}\left(x,t\right)\\\
v_{2}\left(x,t\right)\end{array}\right)=\left(\begin{array}[]{cc}-i\xi&q\left(x,t\right)\\\
\overline{q\left(x,t\right)}&i\xi\end{array}\right)\left(\begin{array}[]{c}v_{1}\left(x,t\right)\\\
v_{2}\left(x,t\right)\end{array}\right),$ (44)

with $\xi\in\mathbb{R}$ being a spectral parameter, is a linear problem
pertinent to the integration of the defocusing cubic nonlinear Schrödinger
(NLS) equation

$i\partial_{t}q\left(x,t\right)=\partial_{x}^{2}q\left(x,t\right)-2\left|q\left(x,t\right)\right|^{2}q\left(x,t\right)$

subject to the initial data $q\left(x,0\right)=q_{0}\left(x\right)$. We refer
to [6] for more details on this matter.

We observe that by setting

$W\left(x,t\right):=\left(\begin{array}[]{cc}e^{i\xi x}&0\\\ 0&e^{-i\xi
x}\end{array}\right)\left(\begin{array}[]{c}v_{1}\left(x,t\right)\\\
v_{2}\left(x,t\right)\end{array}\right),$

Zakharov-Shabat system (44) immediately reduces to

$\partial_{x}W\left(x,t\right)=\left(\begin{array}[]{cc}0&q\left(x\right)e^{2i\xi
x}\\\ \overline{q\left(x\right)}e^{-2i\xi
x}&0\end{array}\right)W\left(x,t\right),$

which is a system of the form (1).

###  3.4 Kubelka-Munk equations

Kubelka-Munk equations is a simple phenomenological model for computing
reflection and transmission optical fluxes without solution of significantly
more complicated radiative transfer equations [9]. Due to their simplicity,
Kubelka-Munk equations have been popular in practice (in paper paint
visibility, see e.g. [3]), they have been extensively studied from modelling
viewpoint and several generalisations have been proposed [13, 14, 15, 16].

We consider the following model equations

$\frac{\textrm{d}}{\textrm{d}x}\left(\begin{array}[]{c}F_{+}\left(x\right)\\\
F_{-}\left(x\right)\end{array}\right)=\left(\begin{array}[]{cc}-K\left(x\right)-S\left(x\right)&S\left(x\right)\\\
-S\left(x\right)&K\left(x\right)+S\left(x\right)\end{array}\right)\left(\begin{array}[]{c}F_{+}\left(x\right)\\\
F_{-}\left(x\right)\end{array}\right),$ (45)

where $F_{+}$, $F_{-}$ are fluxes in positive and negative directions, and $K$
and $S$ are related to absorption and scattering, respectively. Note that,
unlike in the classical model, we take here $K$, $S$ to be dependent on the
optical depth $x$ rather than simply being constants. This generalisation is
expected to be useful since constant scattering and absorption coefficients
are known to be a considerable limitation of the Kubelka-Munk model [3, Sect.
4.5].

The procedure of reduction of (45) to (1) is similar to that performed in
Subsection 3.1. Therefore, we shall omit any detailed calculation.

Let us write

$A\left(x\right):=\left(\begin{array}[]{cc}-K\left(x\right)-S\left(x\right)&S\left(x\right)\\\
-S\left(x\right)&K\left(x\right)+S\left(x\right)\end{array}\right)=S\left(x\right)\left(\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right)-\left[K\left(x\right)+S\left(x\right)\right]\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right),$

and compute

$B\left(x\right):=PA\left(x\right)P^{-1}=S\left(x\right)\left(\begin{array}[]{cc}i&0\\\
0&-i\end{array}\right)-\left[K\left(x\right)+S\left(x\right)\right]\left(\begin{array}[]{cc}0&i\\\
-i&0\end{array}\right)$

with $P$ defined as in (36). Multiplying the both sides of (45) by $P$ and
introducing

$V\left(x\right):=P\left(\begin{array}[]{c}F_{+}\left(x\right)\\\
F_{-}\left(x\right)\end{array}\right),$

we obtain

$V^{\prime}\left(x\right)=B\left(x\right)V\left(x\right).$ (46)

Furthermore, setting

$W\left(x\right):=\left(\begin{array}[]{cc}\exp\left(-i\int_{0}^{x}S\left(\tau\right)\textrm{d}\tau\right)&0\\\
0&\exp\left(i\int_{0}^{x}S\left(\tau\right)\textrm{d}\tau\right)\end{array}\right)V\left(x\right),$
$C\left(x\right):=\left(\begin{array}[]{cc}0&c_{2}\left(x\right)\\\
\overline{c_{2}\left(x\right)}&0\end{array}\right),\hskip
10.00002ptc_{2}\left(x\right):=-i\left[K\left(x\right)+S\left(x\right)\right]\exp\left(-2i\int_{0}^{x}S\left(\tau\right)\textrm{d}\tau\right),$

we arrive at

$W^{\prime}\left(x\right)=C\left(x\right)W\left(x\right),$ (47)

which is a system of the form (1).

##  4 Conclusion

We have introduced a new scalar differential equation of the first order which
is curious for two principal reasons. First, it is, in some sense, the
simplest nonlinear ODE (either with or without a non-homogeneous term), with
the nonlinearity being merely the complex conjugation. Second, this equation
emerges, after appropriate reduction steps, in rather different physical
contexts. Certainly, much more application areas can be identified (e.g.
telegrapher’s equations or Goldstein-Taylor model [5]), but already the
context of the linear Schrödinger equation alone is a good enough motivation
to further study the antilinear ODE
$u^{\prime}\left(x\right)=f\left(x\right)\overline{u\left(x\right)}$. For
example, reduction of matrix-vector manipulations to those involving scalar
quantities already provides a simplification in tedious constructions of
asymptotic-numerical methods, cf. [2, 8]. Therefore, this new reformulation
yields concrete practical advantages. We believe that theoretical aspects of
the mentioned models could benefit from it, too. This might be achievable, for
instance, through newly produced forms of the Prüfer transformation (which is
typically used for studying Sturm-Liouville problems, see e.g. [11, Sect.
5.2]). Furthermore, it is important to identify classes of functions $f$ for
which the antilinear ODE can be solved in a closed form. Here, the Kubelka-
Munk model context hints on the elementary exponential class (note that system
(45) with constant $K$ and $S$ is solvable explicitly). This can be
generalised further since the form of the antilinear ODE is amenable to a
treatment by integral transform methods (unlike other nonlinearities)
typically compatible with an exponential function and combinations thereof.
Finally, the form of the antilinear ODE calls for study of the possible
connection with $d$-bar problems, see e.g. [7]. In this case, an appropriate
extension of the equation to the complex plane may yield a formulation that
eventually produces a closed-form solution due to numerous constructive
results on Hilbert and Riemann-Hilbert problems.

## Acknowledgements

The author is grateful to Austrian Science Fund (FWF) for the support through
the bi-national FWF-project I3538-N32 used for his employment at Vienna
University of Technology (TU Wien). The final version of the manuscript has
benefited from valuable comments of Juliette Leblond to whom the author is
also thankful.

## References

  * [1] M. J. Ablowitz, H. Segur, Solitons and the inverse scattering transform, Society for Industrial and Applied Mathematics (SIAM), 1981.
  * [2] A. Arnold, N. Ben-Abdallah, C. Negulescu, WKB-based schemes for the Schrödinger equation in the semi-classical limit, SIAM J. Math. Anal._,_ 49 (4), 1436–1460, 2011.
  * [3] A. K. R. Choudhury, Principles of Colour and Appearance Measurement - Volume 2 - Visual Measurement of Colour, Colour Comparison and Management, Elsevier Science, 2014.
  * [4] M. Christ, A. Kiselev, WKB Asymptotic Behavior of Almost All Generalized Eigenfunctions for One-Dimensional Schrödinger Operators with Slowly Decaying Potentials, J. Funct. Anal., 179, 426–447, 2001\.
  * [5] H. Dietert, J. Evans, Finding the jump rate for fastest decay in the Goldstein-Taylor model, arXiv:2103.10064, 2021.
  * [6] B. Grébert, T. Kappeler, The defocusing NLS equation and its normal form, European Mathematical Society (EMS), 2014.
  * [7] K. Knudsen, J. Mueller, S. Siltanen, Numerical solution method for the dbar-equation in the plane, J. Comp. Phys., 198 (2), 500–517, 2004.
  * [8] J. Körner, A. Arnold, K. Döpfner, WKB-based scheme with adaptive step size control for the Schrödinger equation in the highly oscillatory regime, J. Comp. Appl. Math., 404, 2022.
  * [9] P. Kubelka, F. Munk, An Article on Optics of Paint Layers (in German; English version available at: http://www.graphics.cornell.edu/~westin/pubs/kubelka.pdf), Z. Tech. Phys, 12, 593–601, 1931.
  * [10] K. Lorenz, T. Jahnke, C. Lubich, Adiabatic integrators for highly oscillatory second-order linear differential equations with time-varying eigendecomposition, BIT _,_ 45, 91–115, 2005.
  * [11] J. D. Pryce, Numerical solution of Sturm-Liouville problems, Oxford University Press, New York, 1993.
  * [12] C. Remling, Spectral theory of canonical systems, De Gruyter, Berlin/Boston, 2018.
  * [13] C. Sandoval, A.D. Kim, Deriving Kubelka-Munk theory from radiative transport, J. Optic. Soc. Am. A, 31 (3), 628–636, 2014\.
  * [14] L. Yang, B. Kruse, Revised Kubelka-Munk theory. I: Theory and application, J. Opt. Soc. Am. A 21 (10), 1933–1941, 2004.
  * [15] L. Yang, B. Kruse, S. J. Miklavcic, Revised Kubelka-Munk theory. II: Unified framework for homogeneous and inhomogeneous optical media, J. Opt. Soc. Am. A 21 (10), 1942–1952, 2004.
  * [16] L. Yang, S. J. Miklavcic, Revised Kubelka-Munk theory. III: A general theory of light propagation in scattering and absorptive media, J. Opt. Soc. Am. A 22 (9), 1866–1873, 2005.


††thanks: These authors have contributed equally to this work.††thanks: These
authors have contributed equally to this work.

# Experimental Determination of Multi-Qubit Ground State via a Cluster Mean-
Field Algorithm

Ze Zhan  Zhejiang Province Key Laboratory of Quantum Technology and Device,
Department of Physics, Zhejiang University, Hangzhou, 310027, China  Chongxin
Run  Zhejiang Province Key Laboratory of Quantum Technology and Device,
Department of Physics, Zhejiang University, Hangzhou, 310027, China  Zhiwen
Zong  Zhejiang Province Key Laboratory of Quantum Technology and Device,
Department of Physics, Zhejiang University, Hangzhou, 310027, China  Liang
Xiang  Zhejiang Province Key Laboratory of Quantum Technology and Device,
Department of Physics, Zhejiang University, Hangzhou, 310027, China  Ying Fei
Zhejiang Province Key Laboratory of Quantum Technology and Device, Department
of Physics, Zhejiang University, Hangzhou, 310027, China  Wenyan Jin  Zhejiang
Province Key Laboratory of Quantum Technology and Device, Department of
Physics, Zhejiang University, Hangzhou, 310027, China  Zhilong Jia  Key
Laboratory of Quantum Information, University of Science and Technology of
China, Hefei, 230026, China  Peng Duan  Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei, 230026, China  Jianlan
Wu<EMAIL_ADDRESS>Zhejiang Province Key Laboratory of Quantum
Technology and Device, Department of Physics, Zhejiang University, Hangzhou,
310027, China  Yi Yin<EMAIL_ADDRESS>Zhejiang Province Key Laboratory of
Quantum Technology and Device, Department of Physics, Zhejiang University,
Hangzhou, 310027, China  Guoping Guo<EMAIL_ADDRESS>Key Laboratory of
Quantum Information, University of Science and Technology of China, Hefei,
230026, China  Origin Quantum Computing, Hefei, 230026, China

###### Abstract

A quantum eigensolver is designed under a multi-layer cluster mean-field (CMF)
algorithm by partitioning a quantum system into spatially-separated clusters.
For each cluster, a reduced Hamiltonian is obtained after a partial average
over its environment cluster. The products of eigenstates from different
clusters construct a compressed Hilbert space, in which an effective
Hamiltonian is diagonalized to determine certain eigenstates of the whole
Hamiltonian. The CMF method is numerically verified in multi-spin chains and
experimentally studied in a fully-connected three-spin network, both yielding
an excellent prediction of their ground states.

Introduction. — At the dawn of a quantum computing era, applications on
quantum simulation and beyond have attracted much attention of the whole
quantum community. For example, mixed quantum-classical algorithms have been
proposed in the goal of solving unaffordable quantum chemistry problems with
quantum computers [1, 2, 3, 5, 4, 6, 7, 8, 9, 10]. A variational quantum
eigensolver (VQE) was successfully implemented in the determination of
electronic states for a hydrogen molecule and multi-atom hydrogen chains [1,
2, 3, 4, 5]. The adiabaticity and shortcut-to-adiabaticity (STA) in analog and
digitized designs [6, 7, 8, 9] can also be used in the quantum eigensolver,
where an eigenstate of the target Hamiltonian is obtained by dragging an
eigenstate of an initial Hamiltonian through an adiabatic or STA trajectory.
Recently, we proposed a ‘leap-frog’ algorithm via the digitized STA and
adiabaticity [10]. Through a segmented trajectory of travelling intermediate
states, our leap-frog method allows an efficient and relibale quantum
eigensolver, as illustrated by our experimental study in H2 and numerical
calculation of hydrogen chains.

In the architecture of quantum computing, the eigenstructure of a
$2^{N}$-dimensional ($2^{N}$-D) Hilbert space can be determined in an
$N$-qubit quantum device. However, the number of quantum gates in a digital
quantum algorithm quickly increases with the number of qubits [12, 11]. In
addition, a multi-qubit quantum gate is realized through a combination of
single- and two-qubit gates [13] but the number of the combining gates
increases with the gate size. The cost of quantum computing increases in
company with the decrease of the fidelity so that a practical quantum
eigensolver is still limited by the system size.

In the fields of physics and chemistry, cluster-based methods have been
applied on various problems [14, 15, 19, 20, 16, 17, 18]. For example, the
concept of block spins was proposed to understand critical phenomena of the
Ising model [14]. In the renormalization group (RG) theory, the critical
exponents are extracted from the scale invariance around a fixed point [15].
The clustering methods are also utilized in the quantum chemistry computation
[19, 20, 16, 17, 18]. In the block renormalization group (BRG) method, the
total Hamiltonian is reconstructed in a compressed Hilbert space built by a
few low-energy block states [16]. In the density matrix renormalization group
(DMRG), the compression of the Hilbert space is realized by the
diagonalization of reduced density matrices [17, 18].

In this paper, we will apply a multi-layer cluster mean-field (CMF) theory
[20] to build a new quantum eigensolver, from which the eigenstructure of a
large-scale system can be reliably and efficiently determined in a much
smaller-scale quantum device. The product states combined from the eigenstates
of reduced cluster Hamiltonians define a compressed Hilbert space, in which
the effective Hamiltonian is diagonalized for the eigensolver. This CMF method
is numerically verified in $N$-spin chains and experimentally implemented in a
fully-connected three-spin system, both yielding high fidelities for the
extracted ground states.

Theory. — In a general multi-electron system, the second quantized Hamiltonian
can be transformed into a multi-spin form,

$\displaystyle
H=g^{(0)}+\sum_{i=1}^{N}\sum_{a=1}^{3}g^{(1)}_{i;a}\sigma_{i}^{a}+\sum_{i,j=1}^{N}\sum_{a,b=1}^{3}g^{(2)}_{ij;ab}\sigma_{i}^{a}\sigma_{j}^{b}+\sum_{i,j,k=1}^{N}\sum_{a,b,c=1}^{3}g^{(3)}_{ijk;abc}\sigma_{i}^{a}\sigma_{j}^{b}\sigma_{k}^{c}+\cdots,$
(1)

through a fermion-to-spin mapping method such as the Bravyi-Kitaev
transformation [21]. Here $\\{\sigma_{i}^{a}\\!=\\!X_{i},Y_{i},Z_{i}\\}$ is
the set of the Pauli matrices acting on spin $i$ and the coefficients
$\\{g^{(0)},g^{(1)},g^{(2)},\cdots\\}$ describe the strengths of (multi)-spin
interactions. To keep its generality, Eq. (1) is allowed to include an
arbitrary $N^{\prime}(\leq N)$-spin interaction.

Figure 1: (a) A schematic diagram of a multi-layer CMF algorithm. (b-c) The
numerical calculation of this CMF method for the $N$-spin systems [the
Hamiltonian in Eq. (3)]: (b) the fidelity $\mathcal{F}^{\mathrm{theo}}_{g}$ of
the ground state and (c) its eigenenergy $E_{g}$. The red circles denote the
numerical results of the CMF while the solid lines denote the exact results.

To extract the exact eigenstates and eigenenergies ($|\Psi_{n}\rangle$ and
$E_{n}$) of the Hamiltonian in Eq. (1), we require a diagonalization tool in a
$2^{N}$-D Hilbert space. Instead, a CMF method can realize an approximate but
reliable eigensolver in a highly compressed space. For simplicity, we assume
that an $N$-spin network is divided into two clusters, each with $N_{A}$ and
$N_{B}(=\\!N\\!-\\!N_{A})$ spins. For a given $N_{A}$, the total choices of
cluster partitioning are $M_{\mathrm{max}}=C^{N}_{N_{A}}$ but a practical
number $M$ can be much smaller than $M_{\mathrm{max}}$.

Next we interpret our CMF method as follows [see Fig. 1(a)]. For a given
$A$-$B$ partition, we inspect the two clusters separately. For cluster $A$,
the rest part of the spin network (cluster $B$) is viewed as its environment.
After a partial trace over a specific $B$-state
$|\varphi^{\alpha}_{B}\rangle$, a reduced Hamiltonian
$H^{\alpha}_{A}=\langle\varphi^{\alpha}_{B}|H|\varphi^{\alpha}_{B}\rangle$ is
constructed and its diagonalization leads to a set of eigenstates
$\\{|\varphi^{i_{\alpha}}_{A}\rangle\\}$ and eigenenergies
$\\{\varepsilon^{i_{\alpha}}_{A}\\}$, i.e.,
$H^{\alpha}_{A}=\sum_{i}\varepsilon^{i_{\alpha}}_{A}|\varphi^{i_{\alpha}}_{A}\rangle\langle\varphi^{i_{\alpha}}_{A}|$.
The same approach can be applied vice versa. With respect to an $A$-state
$|\varphi^{\beta}_{A}\rangle$, the reduced $B$-Hamiltonian
$H^{\beta}_{B}=\langle\varphi^{\beta}_{A}|H|\varphi^{\beta}_{A}\rangle$ is
diagonalized into
$H^{\beta}_{B}=\sum_{j}\varepsilon^{j_{\beta}}_{B}|\varphi^{j_{\beta}}_{B}\rangle\langle\varphi^{j_{\beta}}_{B}|$.
The two sets of product states,
$\\{|\varphi^{i_{\alpha}}_{A}\rangle\otimes|\varphi^{\alpha}_{B}\rangle\\}$
and $\\{|\varphi^{\beta}_{A}\rangle\otimes|\varphi^{j_{\beta}}_{B}\rangle\\}$,
from all the necessary cluster partitions are mixed together to form a basis
set of $\\{|\psi_{\gamma}\rangle\\}$ for a compressed Hilbert space. To
capture a mean-field spirit, we expect that all the states are self-
consistently determined, i.e.,
$\\{|\varphi^{i_{\alpha}}_{A}\rangle\\}=\\{|\varphi^{\beta}_{A}\rangle\\}$ and
$\\{|\varphi^{j_{\beta}}_{B}\rangle\\}=\\{|\varphi^{\alpha}_{B}\rangle\\}$.
Although a regular recursive iteration is divergent if more than one states
are considered, the number of relevant states is in general unchanged. In
practice, we take a limited number of iteration steps. At the final step,
irrelevant states are discarded and the Schmidt orthogonalizations [22] is
used to extract an orthonormal basis set
$\\{|\psi^{\mathrm{S}}_{\gamma}\rangle\\}$. An effective Hamiltonian,

$\displaystyle
H_{\mathrm{eff}}=\sum_{\gamma\gamma^{\prime}}H_{\gamma\gamma^{\prime}}|\psi^{\mathrm{S}}_{\gamma}\rangle\langle\psi^{\mathrm{S}}_{\gamma^{\prime}}|$
(2)

with
$H_{\gamma\gamma^{\prime}}=\langle\psi^{\mathrm{S}}_{\gamma}|H|\psi^{\mathrm{S}}_{\gamma^{\prime}}\rangle$,
is thus defined. The digonalization of $H_{\mathrm{eff}}$ provides a good
estimation of certain eigenstates $|\Psi_{n}\rangle$ and eigenenergies
$E_{n}$. If the number of the product states associated with each cluster
partition is $J$, the dimensionality of the compressed space is $MJ$, which
can be significantly smaller than $2^{N}$. The partition can be subsequently
applied to clusters $A$ and $B$, e.g., $A=a_{1}\oplus
a_{2}=a^{\prime}_{1}\oplus a^{\prime}_{2}=\cdots$, which eventually leads to a
multi-layer CMF algorithm [see Fig. 1(a)]. Relatively speaking, our CMF method
takes a top-down strategy by partitioning a large system into small clusters
while the DMRG takes a bottom-up strategy by extending the system size with
the increment of boundary spins.

Numerical study. — To demonstrate the applicability of this CMF method, we
numerically calculate the ground state $|\Psi_{g}\rangle$ and its eigenenergy
$E_{g}$ of an $N$-spin chain whose Hamiltonian reads [5]

$\displaystyle H=\sum_{i=1}^{N}g_{1}Z_{i}+\sum_{i=1}^{N-1}g_{2}X_{i}X_{i+1}.$
(3)

In our numerical calculation, the chain length is set to be $3\leq N\leq 8$
while the two parameters are fixed at $g_{2}/g_{1}=2$. For each $N$-spin
chain, we only consider two choices of cluster partitioning,
$\\{A=\\{s_{1},s_{2}\\},B=\\{s_{3},\cdots s_{N}\\}\\}$ and
$\\{A^{\prime}=\\{s_{1},\cdots,s_{N-2}\\},B^{\prime}=\\{s_{N-1},s_{N}\\}\\}$,
where $s_{i}$ denotes the $i$-th spin. Taking the first cluster partition as
an example, we show the numerical approach in detail. In the first stage, an
initial $B$-state,
$|\varphi_{B}\rangle\propto\prod\limits_{n=3}^{N}(|+\rangle+|-\rangle)_{n}$,
is used to obtain a reduced $A$-Hamiltonian,

$\displaystyle
H_{A}=\bar{\varepsilon}_{A}+g_{1}Z_{1}+g_{1}Z_{2}+g^{A}_{1}X_{2}+g_{2}X_{1}X_{2}$
(4)

with
$\bar{\varepsilon}_{A}=\sum_{n=3}^{N}g_{1}\langle\varphi_{B}|Z_{n}|\varphi_{B}\rangle$
and $g^{A}_{1}=g_{2}\langle\varphi_{B}|X_{3}|\varphi_{B}\rangle$. Due to the
final goal of calculating $|\Psi_{g}\rangle$, we only select two
$A$-eigenstates, the ground and first excited states of $H_{A}$, i.e.,
$\\{|\varphi_{A}^{\beta}\rangle=|\varphi_{A}^{g}\rangle,|\varphi_{A}^{e}\rangle\\}$.
In the second stage, two $B$-Hamiltonians,

$\displaystyle
H_{B}^{\beta=g,e}=\bar{\varepsilon}_{B}+g^{B}_{1}X_{3}+\sum_{n=3}^{N}g_{1}Z_{n}+\sum_{n=3}^{N-1}g_{2}X_{n}X_{n+1},$
(5)

are extracted with respect to these two $A$-eigenstates. The two parameters
are given by
$\bar{\varepsilon}_{B}=g_{1}\langle\varphi^{\beta}_{A}|Z_{1}+Z_{2}|\varphi^{\beta}_{A}\rangle$
and
$g^{B}_{1}=g_{2}\langle\varphi_{A}^{\beta}|X_{2}|\varphi_{A}^{\beta}\rangle$.
The diagonalization of $H^{\beta=g,e}_{B}$ leads to four (ground and first
excited) $B$-eigenstates, $|\varphi_{B}^{j_{\beta}}\rangle$ with
$j=\\{g,e\\}$. In the third stage, we use these four $B$-states
$|\varphi^{\alpha}_{B}\rangle$ ($\alpha=g_{g},g_{e},e_{g},e_{e}$) as the
environment states and calculate eight $A$-eigenstates
$|\varphi^{i_{\alpha}}_{A}\rangle$ ($i=g,e$). To avoid the divergence of this
recursion, we stop at this stage and discard four crossing terms. The four
remaining products are
$|\varphi_{A}^{i_{\alpha}}\rangle\otimes|\varphi_{B}^{\alpha}\rangle$ with
$\\{i=g,\alpha=g_{g},g_{e}\\}$ and $\\{i=e,\alpha=e_{g},e_{e}\\}$. After
including the four relevant product states from the second cluster partition
$\\{A^{\prime},B^{\prime}\\}$, the Schmidt orthogonalization is applied to
construct an 8-D Hilbert space. The ground state $|\Psi_{g}\rangle$ and its
eigenenergy $E_{g}$ are then determined by the diagonalization of
$H_{\mathrm{eff}}$ in Eq. (2). In the case of $N(>4)$-spin chains, the multi-
layer approach is utilized so that all the calculations are restricted in the
$N^{\prime}(\leq 4)$-spin Hamiltonians. For example, the total 50 two-spin, 5
three-spin and 16 four-spin Hamiltonians are involved for the 8-spin chain.

The numerical results are presented in Figs. 1(b) and 1(c). Here we introduce
a fidelity function,
$\mathcal{F}^{\mathrm{theo}}_{g}\\!=\\!|\langle\Psi^{\mathrm{exact}}_{g}|\Psi^{\mathrm{theo}}_{g}\rangle|^{2}$,
between the CMF result $|\Psi^{\mathrm{theo}}_{g}\rangle$ and the exact state
$|\Psi^{\mathrm{exact}}_{g}\rangle$. As shown in Fig. 1(b), the CMF
predictions are excellent, satisfying $\mathcal{F}^{\mathrm{theo}}_{g}(3\leq
N\leq 8)>99.4\%$. As shown in Fig. 1(c), the accuracy of the ground state
energy $E_{g}$ is even higher ($>99.9\%$) and a good linear dependence is
observed between $E_{g}$ and $N$. As a comparison, we test a DMRG-based
quantum eigensolver and the final state fidelities are around $99.0\%$.

Experimental study. — Next we use a two-qubit device to extract
$|\Psi_{g}\rangle$ and $E_{g}$ of a fully-connected three-spin system as an
experimental demonstration of the CMF algorithm. Due to the restriction of our
current setup, it is difficult for us to reliably explore larger systems which
will be left in the future. Our quantum device is composed of two
superconducting cross-shaped transmon qubits [23, 24, 25]. The ground and
excited states of each qubit are one-to-one mapped onto the spin up and down
states, i.e., $|0\rangle\\!\leftrightarrow\\!|+\rangle$ and
$|1\rangle\\!\leftrightarrow\\!|-\rangle$. The operation points of the two
qubits are $\omega_{a}/2\pi\\!=\\!5.46$ GHz and $\omega_{b}/2\pi\\!=\\!4.92$
GHz, while their anharmonicities are
$\Delta_{a}/2\pi\\!\approx\\!\Delta_{b}/2\pi\\!=\\!-250$ MHz. The relaxation
times are $T_{a;1}\\!=\\!16.1$ $\mu$s and $T_{b;1}\\!=\\!26.5$ $\mu$s, and the
pure dephasing times are $T_{a;\phi}\\!=\\!20$ $\mu$s and
$T_{b;\phi}\\!=\\!45$ $\mu$s. The readout fidelities of the ground and excited
states are $\\{F_{a;0}\\!=\\!99\%,F_{a;1}\\!=\\!93\%\\}$ and
$\\{F_{b;0}\\!=\\!96\%,F_{b;1}\\!=\\!94\%\\}$.

Figure 2: A 4-stage CMF experiment to determine the ground state
$|\Psi_{g}\rangle$ for the 3-spin Hamiltonian in Eq. (6) with
$g_{2}/g_{1}=1.0$ and $g_{3}/g_{1}=0.1$. (a) The 1st-stage eigenenergy
evolution in a 3-segment leap-frog determination of $|\varphi^{g}_{A}\rangle$
for a reduced $A$-Hamiltonian $H_{A}$. (b) The 2nd-stage eigenenergy
evolutions in the digitized STA determination of $|\varphi^{g}_{B}\rangle$ and
$|\varphi^{e}_{B}\rangle$ for the subsequent $B$-Hamiltonian $H_{B}$. (c-d)
The 3rd-stage eigenenergy evolutions in the leap-frog determinations of (c)
$|\varphi^{g_{g}}_{A}\rangle$ and (d) $|\varphi^{e_{e}}_{A}\rangle$ for
$H^{g}_{A}$ and $H^{e}_{A}$, respectively. The cluster partition is shown in
the inset of each panel. In (a), (c) and (d), each cross labels an
intermediate state in the leap-frog algorithm. (e) The 4th-stage eigenenergy
evolution in the VQE determination of $|\Psi_{g}\rangle$ for the 4-D effective
Hamiltonian $H_{\mathrm{eff}}$. The structure of $H_{\mathrm{eff}}$ is
depicted in the inset. In each panel, the symbols denote the experimental
results while the solid horizontal lines label their exact values.

The Hamiltonian of the three-spin system being studied is

$\displaystyle
H=g_{1}(Z_{1}+Z_{2}+Z_{3})+g_{2}(X_{1}X_{2}+X_{2}X_{3})+g_{3}X_{1}X_{2}X_{3},$
(6)

where the three-spin interaction $X_{1}X_{2}X_{3}$ increases the difficulty of
the eigensolver. In this paper, two sets of experiments are performed to
explore the influences of $g_{2}/g_{1}$ and $g_{3}/g_{1}$ separately. In the
first set, we fix $g_{3}/g_{1}=0.1$ and investigate $|\Psi_{g}\rangle$ and
$E_{g}$ upon the change of $g_{2}/g_{1}$. For simplicity, we only consider
three values, $g_{2}/g_{1}=0.1$, 1.0 and 2.0. To visualize our experimental
procedure, we take $g_{2}/g_{1}=1.0$ as an example and provide the stage-by-
stage results in Fig. 2. (i) We treat spins 1 and 2 as cluster $A$ and spin 3
as cluster $B$. With an initial guess of the $B$-state,
$|\varphi_{B}\rangle=|1\rangle$, a reduced $A$-Hamiltonian is obtained as
$H_{A}=H^{0}_{A}+g^{A}_{1}X_{2}+g^{A}_{2}X_{1}X_{2}$ with
$H^{0}_{A}=\bar{\varepsilon}_{A}+g_{1}Z_{1}+g_{1}Z_{2}$. Here the $B$-averaged
parameters are
$\bar{\varepsilon}_{A}=g_{1}\langle\varphi_{B}|Z_{3}|\varphi_{B}\rangle$,
$g^{A}_{1}=g_{2}\langle\varphi_{B}|X_{3}|\varphi_{B}\rangle$ and
$g^{A}_{2}=g_{2}+g_{3}\langle\varphi_{B}|X_{3}|\varphi_{B}\rangle$. The ground
state of $H_{A}$ is experimentally determined by a leap-frog algorithm via the
digitized STA and adiabaticity [26]. With two varying parameters $\lambda_{1}$
and $\lambda_{2}$, the $A$-Hamiltonian is extended to be

$\displaystyle
H_{A}(\lambda_{1},\lambda_{2})=H^{0}_{A}+\lambda_{1}g^{A}_{1}X_{2}+\lambda_{2}g^{A}_{2}X_{1}X_{2}.$
(7)

As shown in Fig. 2(a), we begin with an initial Hamiltonian
$H^{0}_{A}=H_{A}(\lambda_{1}=0,\lambda_{2}=0)$ and prepare its ground state
$|\varphi^{g}_{A}(H^{0}_{A})\rangle=|11\rangle$. A 4-step digitized STA is
applied to drag this state to the ground state
$|\varphi^{g}_{A}(H^{1}_{A})\rangle$ of an intermediate Hamiltonian
$H^{1}_{A}=H_{A}(\lambda_{1}=0,\lambda_{2}=0.1)$. Subsequently, two digitized
adiabatic processes realize an evolution of
$|\varphi^{g}_{A}(H^{1}_{A})\rangle\rightarrow|\varphi^{g}_{A}(H_{A}(\lambda_{1}=0,\lambda_{2}=0.5))\rangle\rightarrow|\varphi^{g}_{A}(H_{A})\rangle$.
The theoretical prediction of the final state fidelity is
$\mathcal{F}^{\mathrm{theo}}_{g}=99.9\%$ while the experimental determination
is at $\mathcal{F}^{\mathrm{exp}}_{g}=98.6\%$. (ii) In the second stage [see
Fig. 2(b)], we input the previous $A$-state $|\varphi^{g}_{A}\rangle$ and
calculate the $B$-Hamiltonian, $H_{B}=H^{0}_{B}+g^{B}_{1}X_{3}$ with
$H^{0}_{B}=\bar{\varepsilon}_{B}+g_{1}Z_{3}$. Here $\bar{\varepsilon}_{B}$ and
$g^{B}_{1}$ are two $A$-averaged parameters. The ground and excited states,
$|\varphi^{g}_{B}(H_{B})\rangle$ and $|\varphi^{e}_{B}(H_{B})\rangle$, are
experimentally determined via the digitized STA from the two initial states
$|\varphi^{g}_{B}(H^{0}_{B})\rangle=|1\rangle$ and
$|\varphi^{e}_{B}(H^{0}_{B})\rangle=|0\rangle$. The experimental fidelities of
these two $B$-eigenstates are $\mathcal{F}^{\mathrm{exp}}_{g}\approx 99\%$.
(iii) In the third stage, the two $B$-states $|\varphi^{g}_{B}\rangle$ and
$|\varphi^{e}_{B}\rangle$ are used to obtain two $A$-Hamiltonians,
$H^{\alpha=g,e}_{A}=\langle\varphi^{\alpha}_{B}|H|\varphi^{\alpha}_{B}\rangle$,
which are extended to the same form $H_{A}(\lambda_{1},\lambda_{2})$ as in Eq.
(7) but the $B$-averaged parameters are updated. As shown in Figs. 2(c) and
2(d), the leap-frog algorithm is also applied to experimentally determine the
ground state $|\varphi^{g_{g}}_{A}\rangle$ of $H^{g}_{A}$ and the first
excited state $|\varphi^{e_{e}}_{A}\rangle$ of $H^{e}_{A}$. The experimental
fidelities are $\mathcal{F}^{\mathrm{exp}}=99.2\%$ and $96.4\%$ while their
theoretical predictions are both $\mathcal{F}^{\mathrm{theo}}=99.8\%$. (iv)
The above iteration stages lead to two product states,
$\\{|\psi_{\gamma=1,2}\rangle=|\varphi^{g_{g}}_{A}\rangle\otimes|\varphi^{g}_{B}\rangle,|\varphi^{e_{e}}_{A}\rangle\otimes|\varphi^{e}_{B}\rangle\\}$.
Following a symmetry argument, the other two product states
$\\{|\psi_{\gamma=3,4}\rangle\\}$ are obtained for the cluster partition of
$A^{\prime}=\\{\mathrm{spins}~{}2,3\\}$ and
$B^{\prime}=\\{\mathrm{spin}~{}1\\}$. The subsequent Schmidt orthogonalization
gives rise to four orthogonal basis states
$\\{\psi^{\mathrm{S}}_{\gamma=1,\cdots,4}\rangle\\}$ and a 4-D effective
Hamiltonian $H_{\mathrm{eff}}$. As shown in Fig. 2(e), the experimental
determination of $|\Psi^{\mathrm{exp}}_{g}\rangle$ is converged over $70\sim
100$ VQE steps, with a high fidelity $\mathcal{F}^{\mathrm{exp}}_{g}=95.4\%$
as compared to the theoretical prediction
$\mathcal{F}^{\mathrm{theo}}_{g}=99.3\%$.

Figure 3: The CMF determination of (a) the fidelity of the ground state
$|\Psi_{g}\rangle$ and (b) the corresponding eigenenergy $E_{g}$ for the
three-spin Hamiltonian in Eq. (6) with a fixed $g_{3}/g_{1}=0.1$ and a varying
$g_{2}/g_{1}$. The red lines and circles denote the numerical and experimental
results via the CMF method while the black lines denote the exact values. In
the inset of (a), the distributions of the $Z$-moment are shown for the
experimentally determined $|\Psi^{\mathrm{exp}}_{g}\rangle$ with
$g_{2}/g_{1}=0.1$, 1.0 and 2.0.

In our numerical calculation of $N$-spin systems, the total 8 product states
are considered in the construction of the compressed Hilbert space. In the
first stage of our experiment, we only consider the ground state of $H_{A}$ so
that the total 4 product states arisen from its first excited state are
excluded. The inset in Fig. 2(e) shows a schematic diagram of the effective
Hamiltonian, from which we find that $|\Psi_{g}\rangle$ can be obtained from
the 3-D or 2-D spaces with $\mathcal{F}^{\mathrm{theo}}_{g}=99.3\%$ and 99.0%.
Thus, a continued compression over the product states is allowed to further
decrease the cost of a CMF eigensolver.

In Fig. 3, we present the experimental results of
$|\Psi^{\mathrm{exp}}_{g}\rangle$ and $E^{\mathrm{exp}}_{g}$ for a fixed
$g_{3}/g_{1}=0.1$ and a varying $g_{2}/g_{1}$ ($=0.1$, 1.0 and 2.0) based on
the CMF algorithm. As compared to the exact ground state, the theoretical
predictions of the state fidelity is excellent
($\mathcal{F}^{\mathrm{theo}}_{g}>99\%$) while the experimental results are
consistently high, $\mathcal{F}^{\mathrm{exp}}_{g}=97.9\%$, $95.4\%$ and
$95.4\%$ [see Fig. 3(a)]. The same behavior is found for the accuracy of
$E^{\mathrm{exp}}_{g}$ [see Fig. 3(b)]. In a simplified scenario of $g_{3}=0$,
this three-spin system prefers ferromagnetism along the $Z$-direction for
$g_{2}/g_{1}\rightarrow 0$ while anti-ferromagnetism along the $X$-direction
in the opposite limit ($g_{2}/g_{1}\rightarrow\infty$). The ground state thus
experiences a transition from $|\Psi_{g}(g_{2}/g_{1}\rightarrow
0)\rangle=|111\rangle$ to
$|\Psi_{g}(g_{2}/g_{1}\rightarrow\infty)\rangle\\!\propto\\!\prod\limits_{i=1}^{3}(|0\rangle-(-1)^{i}|1\rangle)_{i}$.
Here we introduce the total spin moment $M_{j}=\sum_{i=1}^{N}m_{i;j}$, where
$m_{i;j}$ is the magnetic moment of each $i$-th spin along the
$j(=X,Y,Z)$-direction. In the parameter range in our experiment, the
entanglement of $|\Psi_{g}\rangle$ increases with $g_{2}/g_{1}$, indicated by
a broadening distribution of the $Z$-moment $M_{Z}$ in the inset of Fig. 3(a).

Figure 4: The CMF determination of (a) the fidelity of the ground state
$|\Psi_{g}\rangle$ and (b) the corresponding eigenenergy $E_{g}$ for the
three-spin Hamiltonian in Eq. (6) with a fixed $g_{2}/g_{1}=2.0$ and a varying
$g_{3}/g_{1}$. The red lines and circles denote the numerical and experimental
results via the CMF method while the black lines denote the exact values. In
the inset of (a), the distributions of the $Z$-moment are shown for the
experimentally determined $|\Psi^{\mathrm{exp}}_{g}\rangle$ with
$g_{3}/g_{1}=0.1$, 1.0 and 2.0.

In our second set of experiments, we fix $g_{2}/g_{1}=2.0$ and consider three
values of $g_{3}/g_{1}=0.1$, 1.0 and 2.0. As shown in the inset of Fig. 4(a),
the increase of $g_{3}/g_{1}$ also leads to an extensive distribution of the
$Z$-moment. The same CMF algorithm as in Fig. 2 reliably determines the ground
states. The experimental results of the state fidelities are presented in Fig.
4(a), satisfying $\mathcal{F}^{\mathrm{exp}}_{g}=95.4\%$, $94.0\%$ and
$95.3\%$ for the three input parameters. The accuracy of the experimentally
extracted eigenenergy $E^{\mathrm{exp}}_{g}$ follows the same trend [see Fig.
4(b)].

Summary. — In this paper, we apply a multi-layer CMF method to design a new
quantum eigensolver so that the eigenstates of a large-scale quantum system
can be determined by a series of quantum computations over its clusters. For a
pre-selected cluster, certain eigenstates of its reduced Hamiltonian are
extracted via a quantum algorithm after a partial average over an eigenstate
of the environment cluster. The products of eigenstates from different
clusters are used to construct a compressed Hilbert space, in which the
effective Hamiltonian is digonalized to determine certain eigenstates of the
whole Hamiltonian. This CMF method is numerically verified in the $N$-spin
chains with two-spin interactions. For the condition of $3\leq N\leq 8$, the
CMF calculations in the $2^{N^{\prime}}(N^{\prime}\leq 4)$-D spaces provide an
excellent prediction on the ground state of the $2^{N}$-D Hamiltonian with
$\mathcal{F}^{\mathrm{theo}}_{g}>99.4\%$. This CMF method is further
experimentally studied in the 3-spin chain with both two- and three-spin
interactions. Under various parameter combinations, the experimental
determination of the ground state via the CMF method is consistently high,
satisfying $\mathcal{F}^{\mathrm{exp}}_{g}\gtrsim 95\%$. The studies of the
ground states in this paper can be straightforwardly extended to the excited
states. With a quick size increment of quantum devices [27, 28], the CMF
method shows its promise to sufficiently large-scale Hilbert spaces.

###### Acknowledgements.

The work reported here was supported by the National Key Research and
Development Program of China (Grant No. 2019YFA0308602, No. 2016YFA0301700),
the National Natural Science Foundation of China (Grants No. 12074336, No.
11934010, No. 11775129), the Fundamental Research Funds for the Central
Universities in China (2020XZZX002-01), and the Anhui Initiative in Quantum
Information Technologies (Grant No. AHY080000). Y.Y. acknowledge the funding
support from Tencent Corporation. This work was partially conducted at the
University of Science and Technology of China Center for Micro- and Nanoscale
Research and Fabrication.

## References

  * [1] A. Peruzzo, J. McClean, P. Shadbolt, M. Yung, X. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nat. Commun. 5, 4213 (2014).
  * [2] P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, and N. Ding et al., Phys. Rev. X 6, 031007 (2016).
  * [3] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * [4] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi, Phys. Rev. X. 8, 011021 (2018).
  * [5] Google AI Quantum, Science 369, 1084 (2020).
  * [6] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, D. Preda, Science 292, 472 (2001).
  * [7] R. Barends, A. Shabani, L. Lamata, J. Kelly, A. Mezzacapo, U. L. Heras, R. Babbush, A. G. Fowler, B. Campbell, and Yu Chen et al., Nature 534, 222 (2016).
  * [8] X. Chen, I. Lizuain, A. Ruschhaupt, D. Guéry-Odelin, and J. G. Muga, Phys. Rev. Lett. 105, 123003 (2010).
  * [9] N. N. Hegade, K. Paul, Y. Ding, M. Sanz, F. Albarran-Arriagada, E. Solano, and X. Chen, Phys. Rev. Appl. 15, 024038 (2021).
  * [10] Z. Zhan, C. Run, Z. Zong, L. Xiang, Y. Fei, Z. Sun, Y. Wu, Z. Jia, P. Duan, J. Wu, Y. Yin, and G. Guo, arXiv: 2103.06098.
  * [11] S. Lloyd, Science 273, 1073 (1996).
  * [12] R. Barends, L. Lamata, J. Kelly, L. García-Álvarez, A. G. Fowler, A. Megrant, E Jeffrey, T. C. White, D. Sank, and J. Y. Mutus et al., Nat. Commun. 6, 7654 (2015).
  * [13] S. Lloyd, Phys. Rev. Lett. 75, 346 (1995).
  * [14] L. P. Kadanoff, Physics 2, 263 (1966).
  * [15] K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
  * [16] S. D. Drell, M. Weinstein, and S. Yankielowicz, Phys. Rev. D 14, 487 (1976).
  * [17] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * [18] U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
  * [19] D. Yamamoto, Phys. Rev. B 79, 144427 (2009).
  * [20] C. A. Jiménez-Hoyos and G. E. Scuseria, Phys. Rev. B 92, 085101 (2015).
  * [21] J. T. Seeley, M. J. Richard, and P. J. Love, J. Chem. Phys. 137, 224109 (2012).
  * [22] Å. Björck, BIT 7, 1 (1967).
  * [23] A. A. Houck, D. I. Schuster, J. M. Gambetta, J. A. Schreier, B. R. Johnson, J. M. Chow, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin and R. J. Schoelkopf, Nature 449, 328 (2007).
  * [24] R. Barends, J. Kelly, A. Megrant, D. Sank, E. Jeffrey, Y. Chen, Y. Yin, B. Chiaro, J. Mutus, and C. Neill et al., Phys. Rev. Lett. 111, 080502 (2013).
  * [25] L. Xiang, Z. Zong, Z. Sun, Z. Zhan, Y. Fei, Z. Dong, C. Run, Z. Jia, P. Duan, J. Wu, Y. Yin, and G. Guo, Phys. Rev. Appl. 14, 014099 (2020).
  * [26] See Supplementary Material at [URL will be inserted by publisher] for a discussion of the leap-frog algorithm.
  * [27] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, and D. A. Buell et al., Nature 574, 505 (2019).
  * [28] M. Gong, S. Wang, C. Zha, M. Chen, H. Huan, Y. Wu, Q. Zhu, Y. Zhao, S. Li, and S. Guo et al., Science 372, 948 (2021).


# All-electrical detection of skyrmion lattice state and chiral surface twists

A. Aqeel  Department of Physics, Technical University Munich, 85748 Garching,
Germany  M. Azhar  Zernike Institute for Advanced Materials, University of
Groningen, Nijenborgh 4, 9747 AG Groningen, The Netherlands  N. Vlietstra
Department of Physics, Technical University Munich, 85748 Garching, Germany
Walther-Mei$\ss$ner-Institute, 85748 Garching, Germany  A. Pozzi  Zernike
Institute for Advanced Materials, University of Groningen, Nijenborgh 4, 9747
AG Groningen, The Netherlands  J. Sahliger  Department of Physics, Technical
University Munich, 85748 Garching, Germany  H. Huebl  Walther-Mei$\ss$ner-
Institute, 85748 Garching, Germany  Department of Physics, Technical
University Munich, 85748 Garching, Germany  Munich Center for Quantum Science
and Technology (MCQST), Schellingstr. 4, D-80799 Munich,Germany  T. T. M.
Palstra  Zernike Institute for Advanced Materials, University of Groningen,
Nijenborgh 4, 9747 AG Groningen, The Netherlands  C. H. Back  Department of
Physics, Technical University Munich, 85748 Garching, Germany  Munich Center
for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799
Munich,Germany  M. Mostovoy  Zernike Institute for Advanced Materials,
University of Groningen, Nijenborgh 4, 9747 AG Groningen, The Netherlands

(August 27, 2024)

###### Abstract

We study the high-temperature phase diagram of the chiral magnetic insulator
Cu2OSeO3 by measuring the spin-Hall magnetoresistance (SMR) in a thin Pt
electrode. We find distinct changes in the phase and amplitude of the SMR
signal at critical lines separating different magnetic phases of bulk
Cu2OSeO3. The skyrmion lattice state appears as a strong dip in the SMR phase.
A strong enhancement of the SMR amplitude is observed in the conical spiral
state, which we explain by an additional symmetry-allowed contribution to the
SMR present in non-collinear magnets. We demonstrate that the SMR can be used
as an all-electrical probe of chiral surface twists and skyrmions in magnetic
insulators.

Magnetic skyrmions are nano-scale spin-swirling objects of great interest due
to their small size and topological protection Nagaosa and Tokura (2013); Fert
_et al._ (2017). These magnetic solitons are being considered as promising
candidates for information bits in ultra compact thin film memory devices
where surfaces and interfaces play a crucial role. Recent theoretical studies
of chiral magnets showed that their surface magnetization can be strongly
modified by the bulk Dzyaloshinskii-Moriya interactions (DMI), leading to a
surface structure which deviates from the bulk skyrmion lattice. These so-
called chiral surface twists are expected to be prominent in the field-
polarized magnetic state restraining the magnetization to fully saturate near
the surface along the applied field direction Rybakov  _et al._ (2013). They
also affect helicity and stability of three-dimensional skyrmions Rybakov  _et
al._ (2013); Meynell  _et al._ (2014); Rybakov  _et al._ (2015); Leonov  _et
al._ (2016) and first experimental evidence for the existence of these chiral
surface twists has been recently reported for skyrmion tubes aligned
perpendicular to the sample surface Zhang  _et al._ (2018). Chiral twists were
also observed in thin film multilayers Legrand  _et al._ (2018); Zhang  _et
al._ (2020) potentially relevant for skyrmionics applications.

In this paper, we report all-electrical detection of skyrmions in Cu2OSeO3
(CSO). Close to spin-ordering temperature, this chiral magnetic insulator with
P$2_{1}3$ symmetry hosts crystals of three-dimensional tubes formed by stacked
Bloch-type skyrmions Seki  _et al._ (2012). The effect of the chiral twist on
the electrical fingerprint of the skyrmion tubes, which in our device are
parallel to the CSO surface, has not been studied so far.

The detection of these surface modifications is difficult, yet important for
integration of skyrmions in multi-layered devices. Here, we systematically
track the surface magnetization by a detailed investigation of the angular
dependence of the spin Hall magnetoresistance (SMR) in the high-temperature
part $(40\leavevmode\nobreak\ {\rm K}-60\leavevmode\nobreak\ {\rm K})$ of the
CSO phase diagram.

Figure 1: (a) Schematic representation of the device configuration used to
measure the spin-Hall magnetoresistance in the Pt$|$CSO bilayer system, with
an applied magnetic field $H$ at an angle $\alpha$ to $x$, the direction of
the current. The skyrmion tubes are represented by color-coding the
magnetization projection along the applied magnetic field direction (red to
blue). Top (b) and cross-sectional (c) view of the surface twists in the
field-polarized state. The increase in the twist angle between the
magnetization and the applied field near the surface is indicated by color
(red to green).

SMR has been used in the past to electrically detect the in-plane
magnetization of collinear magnetic insulators Vlietstra  _et al._ (2013);
Nakayama  _et al._ (2013); Althammer  _et al._ (2013), non-collinear magnets
Aqeel  _et al._ (2015); Ganzhorn  _et al._ (2016); Aqeel  _et al._ (2016,
2017), and antifermagnetic materials Ji  _et al._ (2017); Hoogeboom  _et al._
(2017); Wang  _et al._ (2017); Fischer  _et al._ (2018); Lebrun  _et al._
(2019). In the SMR, the resistance of a heavy metal (Pt) is sensitive to the
magnetization direction of an adjacent magnetic layer. In collinear magnets,
the SMR is theoretically approximated by the magnetic moment density at the
interface Jia  _et al._ (2011); Chen  _et al._ (2013). In addition to a
longitudinal resistance change, a transverse voltage arises given by

$V_{\rm SMR}\propto m_{x}m_{y},$ (1)

where $m_{x}$ and $m_{y}$ are the in-plane components of the unit vector
$\bm{m}$ describing the magnetization direction. For non-collinear magnets
with a spin relaxation length, $\xi$, much smaller than the typical length
scale of variation of $\bm{m}(\bm{x})$, the right-hand side of Eq.(1) is
replaced by its average over the interface Aqeel  _et al._ (2016):

$V_{\rm SMR}\propto\langle m_{x}m_{y}\rangle.$ (2)

When the external magnetic field $H$ (larger than the saturation magnetic
field $H_{\rm c2}$) is rotated in the plane of the normal metal$|$magnet
interface, the SMR voltage, $V_{\rm SMR}$, follows a sinusoidal angular
dependence with a periodicity of 180∘ (see Fig. 2(a)). This behavior of the
SMR is expected when the magnetization is fully aligned with the applied
magnetic field. However, any in-plane tilt of the magnetization away from $H$
will give rise to an additional phase, $\phi$, in the angular dependence of
$V_{\rm SMR}$ making this physical parameter very sensitive to the magnetic
structure at the interface.

Figure 2: (a) Angular dependence of the SMR signal, $V_{\rm SMR}$, measured at
$\mu_{0}H=600$ mT in a planar-Hall geometry shown in Fig. 1(a). $V_{\rm
SMR}^{\rm ampl}$ and $\phi$ are the amplitude and phase of $V_{\rm SMR}$,
respectively. The solid lines show the $V_{\rm SMR}^{\rm
ampl}\sin$(2$(\alpha-\phi)$) fit with $\phi=0^{\circ}$ (black) and
$\phi=90^{\circ}$ (red). (b) and (c) show the temperature dependence of
$V_{\rm SMR}^{\rm ampl}$ and $\phi$, respectively. Here, the shaded area
corresponds to the skyrmion lattice state. (d) shows the contour map of phase
$\phi$ versus magnetic field and temperature combined with the results of
broadband ferromagnetic resonance. The hollow and solid circles represent the
phase boundaries extracted from the magnetic resonance data measured by
applying $\bm{H}$ along the [1 -1 1] and [-2 -1 1] crystallographic directions
of CSO, respectively. The magnetic phase boundaries are highlighted by dashed
lines. We identify helical (H), conical (C) ($H_{\rm c1}$ is the low field
boundary), field-polarized ferrimagnetic state (FP) ($H_{\rm c2}$ is the low
field boundary) and the skyrmion lattice state (SkL) ($H_{\rm A1}$ and $H_{\rm
A2}$ are low and high field boundaries). (e) shows the line scans of $\phi$ as
function of applied field $H$ at two different temperatures.

It is worth noting that previous SMR studies of CSO were performed outside the
stability region of the Skyrmion lattice (SkL) phase Aqeel  _et al._ (2016,
2017). The measurements reported here focus on the SkL. We show that the SMR
sensitively measures the deviation of the local interface magnetization from
$\bm{H}$, allowing us to probe the highly non-collinear SkL spin texture
electrically. The obtained results provide new insights into electrical
detection of skyrmions and their surface deformations in magnetic insulators.

SMR measurements are performed in the device configuration shown in Fig. 1(a),
where a 5 nm thick Pt Hall cross is structured onto an oriented polished CSO
cuboid (dimensions: 4$\times$4$\times$2 mm3). A low frequency ac-current
($f=17$ Hz) is applied through Pt along the x axis. The SMR is detected in the
transverse configuration (along the y axis) by measuring the first harmonic
voltage response $V_{1}$ using a lock-in amplifier Vlietstra  _et al._ (2014)
(see Supplementary Material for more details). The angular dependence of the
SMR is recorded by rotating the applied magnetic field within the xy-plane
parallel to the Pt/CSO interface. Here, the assumption that the SKL follows
the field direction, has been recently confirmed in CSO for an in-plane field
rotation in a different crystal orientation Zhang  _et al._ (2020). A typical
example of the angular dependence of the SMR data measured at 5 K in the
ferrimagnetic state of CSO ($\mu_{0}H=600$ mT) is shown in Fig. 2(a). The
amplitude $V_{\rm SMR}^{\rm ampl}$ and phase $\phi$ are extracted by fitting
the data with $V_{\rm SMR}=V_{\rm SMR}^{\rm ampl}\sin(2(\alpha-\phi))$. Here,
$\alpha$ is the angle between the applied current $I$ and the applied magnetic
field, $\bm{H}$, as defined in Fig. 1(a). The phase $\phi$ is zero, if
magnetization is parallel to the field. The $\sin(2\alpha)$ dependence with
zero phase has been observed in the ferrimagnetic state of CSO at 5K (see Fig.
2(a)). On the other hand, in the helical and low-field conical states, the
magnetic moments align almost perpendicular to the propagation vector of the
spiral, resulting in $\phi\sim 90^{\circ}$ Aqeel  _et al._ (2016).

To identify the skyrmion pocket, we recorded the angular dependence of the SMR
using two sets of measurement protocols: (i) at fixed fields as a function of
temperature (T-scan) and (ii) at fixed temperature, recording the SMR angular
dependence at various magnetic fields strengths (H-scan). Figures 2(b) and
2(c) show the T-scan of $\phi$ and $V_{\rm SMR}^{\rm ampl}$ extracted from the
angular dependence of the SMR signal as exemplified in Fig. 2(a). At
$\mu_{0}H=20$ mT outside the skyrmion pocket, CSO is expected to be in the
conical spiral state, in which $\phi\sim 90^{\circ}$ is expected. Indeed, we
observe that in a wide temperature range $\phi\sim 90^{\circ}$ except near
$T_{\rm c}$ where a clear dip in the phase $\phi$ is found (labeled as SkL in
Fig. 2(b)). A similar signature dip is observed in $V_{\rm SMR}^{\rm ampl}$
(see Fig. 2(c); note the logarithmic voltage scale).

Figure 3: Field dependence of the spin-Hall magnetoresistance (a) amplitude,
$V_{\rm SMR}^{\rm ampl}$, and (b) phase, $\phi$, in the conical spiral (C) and
collinear field-polarized (FP) states calculated for each of the three terms
in Eq. (3) separately: blue ($A_{1}$ and $\phi_{1}$, due to the 1st term),
magenta ($A_{2}$ and $\phi_{2}$, due to the 2nd term) and green lines ($A_{3}$
and $\phi_{3}$, due to the 3rd term). Field dependence of (c) $V_{\rm
SMR}^{\rm ampl}$ and (d) $\phi$ in the C and FP states calculated numerically
using Eq. (3) with $bQ/a=-1.59$ and $c/a=0$ (solid blue lines). The dashed
line shows the effect of the interfacial DMI, for $D_{\rm int}=-0.30J$. The
dotted blue line is obtained by considering only the first term in Eq. (3).
Open circles are the experimental data taken at 51 K. Here, the critical
field, $H_{\rm c2}$ = 56 mT. Field-dependence of (e) $V_{\rm SMR}^{\rm ampl}$
and (f) $\phi$ calculated numerically using Eq. (3) with $bQ/a=-1.84$ and
$c/a=0$ for the C and FP states (blue solid lines) and for the skyrmion
lattice (SkL) state (red solid lines). The dashed blue (red) lines show the
effect of the interfacial DMI of strength $D_{\rm int}=-0.40J$ for the C (SkL)
state. Open circles are results of experimental measurements at 55 K. The
color-bars at the $x$ axes of the figures indicate the magnetic phases,
deduced from FMR measurements for the experimental data: the C, SkL and FP
state is indicated by yellow, red and white color, respectively. The critical
fields determined from FMR measurements at 55 K are: $H_{\rm c2}=40$ mT, and
the boundaries of the SkL, $H_{\rm A1}$ and $H_{\rm A2}$, are 15 mT and 23 mT,
respectively.

To outline the boundaries of the skyrmion pocket, we constructed the high-
temperature phase diagram of CSO by mapping the phase $\phi$ as a function of
field and temperature using H-scans. Figure 2(d) shows the obtained phase
diagram along with the results of ferromagnetic resonance (FMR) measurements
which are used here to identify the phase boundaries of the different magnetic
states of the CSO crystal. To measure the FMR response, a broad band spin-wave
spectroscopy technique Schwarze  _et al._ (2015) is used (see the
Supplementary Material for details). The phase boundaries of the FMR data
(shown by circles in Fig. 2(d)) agree well with those extracted from the SMR
phase $\phi$. $V_{\rm SMR}^{\rm ampl}$ also exhibits distinct anomalies at the
magnetic phase transitions, though less pronounced than those in the SMR phase
$\phi$ (see Fig. S3 in the Supplementary Material). From the experimentally
obtained phase diagram shown in Fig. 2(d), the main observations are as
follows: (i) the SkL state is clearly distinguishable from other magnetic
states of CSO, and (ii) close to $T_{\rm c}$, the phase $\phi$ remains non-
zero above $H_{c2}$. To comprehend these findings, two H-scans of $\phi$ at
two temperatures are shown in Fig. 2(e). At 54 K, the phase $\phi$ is almost
90∘ in the helical state and decreases in the conical state of CSO; the lowest
values of the phase $\phi$ are observed in the FM state. At 55 K, the phase
$\phi$ follows the same trend as observed at 54 K and in addition shows a
significant change in the SkL state of CSO. Importantly, in the collinear
state ($H\geq H_{\rm c2}$), $\phi$ remains non-zero ($\sim 50^{\circ}$ at
$\mu_{0}H$ = 50 mT).

A nonzero SMR phase above $H_{\rm c2}$ provides direct evidence for a chiral
surface twist, i.e. the tilt of the magnetization at the interface away from
the magnetic field direction Rybakov  _et al._ (2013); Meynell  _et al._
(2014); Rybakov  _et al._ (2015); Leonov  _et al._ (2016). The observed field-
dependence of $\phi$ can be understood using Eq. (2) and taking into account
the modification of the conical spiral near the surface. We numerically
calculated this modification by minimizing the energy of a classical spin
model describing CSO with periodic boundary conditions along the in-plane
directions and open boundary conditions along the surface normal direction
(dotted blue line in Fig. 3(d); see Supplementary Material for details of the
calculations). This approach reproduces the experimentally observed decrease
of the SMR phase $\phi$ with increasing magnetic field in the conical spiral
state and a nonzero value of $\phi$ at the transition to the field-induced
collinear state: due to the rotation of spins at the interface away from the
magnetic field direction, $\phi=$$$ at $H=H_{\rm c2}$ and then slowly
decreases with increasing magnetic field for $H>H_{\rm c2}$, in agreement with
our observations (cf. measured data (open circles) with dotted blue line in
Fig. 3(d)).

However, this approach fails to explain the observed field-dependence of the
SMR amplitude. The calculated $V_{\rm SMR}^{\rm ampl}$ reaches its maximal
value at $H_{\rm c2}$ and stays constant at higher fields (dotted blue line in
Fig. 3(c)). In contrast, the SMR amplitude measured above 50 K is relatively
small at $H\geq H_{\rm c2}$. It reaches a maximum at rather low magnetic
fields and then decreases with increasing field (see Fig. 3(c)). The origin of
this discrepancy can be traced back to Eq. (2) that only involves in-plane
magnetization components which become on average larger as the strength of the
magnetic field parallel to the interface increases. As a result, the
calculated SMR amplitude grows, as $H$ tends towards $H_{\rm c2}$. This
discussion shows that the surface twist alone cannot explain the high-
temperature SMR data and motivates us to add new phenomenological terms to Eq.
(2) allowed by symmetry of the polycrystalline Pt film (rotations around the
normal to the film and vertical mirrors):

$\begin{split}\rho_{\rm
SMR}^{\perp}=\quad&a<m_{x}m_{y}>+b<m_{z}\overset{\leftrightarrow}{\partial_{y}}m_{x}>\\\
\\\ &+c<m_{z}^{2}m_{x}m_{y}>,\end{split}$ (3)

where $\rho_{\rm SMR}^{\perp}$ is the transverse SMR resistivity. The first
term coincides with Eq. (2). The second term is the contribution proportional
to the magnetization gradient along the in-plane direction:
$m_{z}\overset{\leftrightarrow}{\partial_{y}}m_{x}=m_{z}(\partial_{y}m_{x})-(\partial_{y}m_{z})m_{x}$.
This term is sensitive to the presence of the helical spin spiral modulation
at the interface, and is proportional to the spiral wavevector,
$Q=2\pi/\lambda$, where $\lambda$ is the period of the spiral. The last term
proportional to the fourth power of magnetization can originate from higher-
order spin torques Hanke  _et al._ (2020) (the higher-order terms in $\bm{m}$
are omitted, for simplicity). The microscopic mechanism behind the second term
in Eq. (3) is unclear. According to SMR theory Chen  _et al._ (2013), it
should be a $\xi Q$ correction to the first term, $\xi\sim 1$ nm being the
relaxation length in Pt, whereas our experiment suggests that the
contributions of first and second terms of Eq. (3) are comparable.

The field-dependence of these three contributions to the SMR amplitude and
phase is shown in Figs. 3(a,b). Both the second and the third terms can
explain the observed decrease of $V_{\rm SMR}^{\rm ampl}$ with increasing
magnetic field. However, the third term gives rise to an extra minimum, which
is not observed in experiment. The best fit of the field-dependence of both
the amplitude and phase of SMR is obtained for $bQ/a=-1.59$, $c=0$ at 51 K and
$bQ/a=-1.84$, $c=0$ at 55 K (solid blue lines in Fig. 3(c-f)). The surprising
conclusion is that the contribution of the new term proportional to the
magnetization gradient (the second term in Eq. (3)) to the SMR of CSO is
comparable to that of the first term in Eq. (3) that works well for collinear
magnets.

Although the additional SMR term improves the overall agreement between the
calculated and experimental SMR curves, the conical spiral state in the
theoretical plots seems to disappear at a field that is lower than $H_{\rm
c2}$ deduced from our SMR and FMR measurements (see solid blue lines in Fig.
3(c-f)). This mismatch reflects an interesting effect found by numerical
simulations: in addition to the rotation of spins around the surface normal,
present in both the collinear and conical states, the conical angle at the
interface is smaller than that in bulk. Moreover, it becomes very small at
fields well below $H_{\rm c2}$, so that, in practice, the interface has a
lower critical field than the bulk (see Figs. S6 and S7 in the Supplementary
Material).

Our simulations also reveal strong sensitivity of the SMR to additional spin
interactions at the interface, such as the interfacial Dzyaloshinskii-Moriya
interaction (iDMI) enabled by inversion symmetry breaking at Pt$|$CSO
interface and the strong spin-orbit coupling of Pt. The iDMI with a positive
(negative) constant, $D_{\rm int}$, decreases (increases) the cone angle of
the spiral at the interface, leaving the surface twist angle practically
unchanged. Hence, it affects the effective critical field at the interface,
which can solve the above mentioned problem. Figures 3(c-f) show that iDMI can
improve the agreement between the theoretically calculated SMR amplitude and
phase (blue dashed line) and the experimental data. The calculations are done
for $D_{\rm int}=-0.40J$ at 55 K and $D_{\rm int}=-0.30J$ at 51 K, where $J$
is the nearest-neighbor Heisenberg exchange constant.

Finally, we discuss the sudden changes in the experimentally measured
amplitude and phase of the SMR associated with the intervening skyrmion
lattice phase. Figures 3(e,f) show the field dependence of $V_{\rm SMR}^{\rm
ampl}$ and $\phi$ in the SkL state calculated numerically using Eq. (3) with
and without the iDMI (dashed and solid red lines, respectively). Although
minimization of the energy of stable and metastable spin configurations at
zero temperature does not allow us to obtain boundaries of the SkL phase near
$T_{c}$, our calculations show that the SMR in the SkL state can be very
different from that in the conical spiral state. In particular, the calculated
phase $\phi$ in the SkL state is smaller than that in the conical spiral state
at low magnetic fields, which agrees with the observed dip in the SMR phase.
In addition, we find a sudden jump in both the phase and the amplitude of the
SMR in the SkL state related to an abrupt change of the skyrmion crystal spin
configuration at the surface (See Fig. S9 in Supplementary Materials). At high
fields, skyrmions are repelled from the edge and spins at the interface are
collinear and nearly parallel to the surface. At low fields, the skyrmion
centers are residing at the interface, which makes their topological charge
smaller than 1 Keesman  _et al._ (2015). The field at which this change occurs
is strongly affected by the interfacial DMI (red dashed lines in Fig. 3).

To summarize, our SMR measurements provide a clear evidence for a surface
twist in the conical and collinear phases near $T_{c}$ and make it possible to
all-electrically outline the boundaries of the skyrmion crystal phase, at
which both the amplitude and the phase of the SMR show a profound
discontinuity. The theoretical description of SMR for non-collinear chiral
magnets, such as CSO, is more involved than that for collinear ferromagnets.
First, SMR is affected by subtle interfacial effects, such as the difference
between the surface and bulk critical fields. Second, a new phenomenological
term in the SMR expression was required to reproduce the field-dependence of
the SMR amplitude. In addition, the SMR was found to be sensitive to
interfacial interactions. Our experimental and theoretical results
conclusively show that the SMR, as a probe of the surface magnetization, is an
important tool for the all-electrical detection of magnetic phase transitions
in chiral magnets.

###### Acknowledgements.

We would like to thank B. J. van Wees, T. Kuschel and J. Bass for useful
discussions. This work has been funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) via project number 107745057 - TRR 80, via
SPP2137, as well as via Germany’s Excellence Strategy – EXC-2111 – 390814868,
Vrije FOM-programma ‘Skyrmionics’ and the Zernike Institute for Advanced
Materials.

## References

  * Nagaosa and Tokura (2013) N. Nagaosa and Y. Tokura, Nature Nanotechnology 8, 899 (2013).
  * Fert  _et al._ (2017) A. Fert, N. Reyren, and V. Cros, Nature Reviews Materials 2, 17031 (2017).
  * Rybakov  _et al._ (2013) F. N. Rybakov, A. B. Borisov, and A. N. Bogdanov, Phys. Rev. B 87, 094424 (2013).
  * Meynell  _et al._ (2014) S. A. Meynell, M. N. Wilson, H. Fritzsche, A. N. Bogdanov, and T. L. Monchesky, Phys. Rev. B 90, 014406 (2014).
  * Rybakov  _et al._ (2015) F. N. Rybakov, A. B. Borisov, S. Blügel, and N. S. Kiselev, Phys. Rev. Lett. 115, 117201 (2015).
  * Leonov  _et al._ (2016) A. O. Leonov _et al._ , Phys. Rev. Lett. 117, 087202 (2016).
  * Zhang  _et al._ (2018) S. L. Zhang, G. van der Laan, W. W. Wang, A. A. Haghighirad, and T. Hesjedal, Phys. Rev. Lett. 120, 227202 (2018).
  * Legrand  _et al._ (2018) W. Legrand _et al._ , Science Advances 4 (2018), 10.1126/sciadv.aat0415.
  * Zhang  _et al._ (2020) S. Zhang _et al._ , Nano Letters 20, 1428 (2020).
  * Seki  _et al._ (2012) S. Seki, X. Z. Yu, S. Ishiwata, and Y. Tokura, Science 336, 198 (2012).
  * Vlietstra  _et al._ (2013) N. Vlietstra, J. Shan, V. Castel, B. J. van Wees, and J. Ben Youssef, Phys. Rev. B 87, 184421 (2013).
  * Nakayama  _et al._ (2013) H. Nakayama _et al._ , Phys. Rev. Lett. 110, 206601 (2013).
  * Althammer  _et al._ (2013) M. Althammer _et al._ , Phys. Rev. B 87, 224401 (2013).
  * Aqeel  _et al._ (2015) A. Aqeel _et al._ , Phys. Rev. B 92, 224410 (2015).
  * Ganzhorn  _et al._ (2016) K. Ganzhorn _et al._ , Phys. Rev. B 94, 094401 (2016).
  * Aqeel  _et al._ (2016) A. Aqeel _et al._ , Phys. Rev. B 94, 134418 (2016).
  * Aqeel  _et al._ (2017) A. Aqeel, M. Mostovoy, B. J. van Wees, and T. T. M. Palstra, Journal of Physics D: Applied Physics 50, 174006 (2017).
  * Ji  _et al._ (2017) Y. Ji _et al._ , Applied Physics Letters 110, 262401 (2017).
  * Hoogeboom  _et al._ (2017) G. R. Hoogeboom, A. Aqeel, T. Kuschel, T. T. M. Palstra, and B. J. van Wees, Applied Physics Letters 111, 052409 (2017).
  * Wang  _et al._ (2017) H. Wang _et al._ , Journal of Applied Physics 122, 083907 (2017).
  * Fischer  _et al._ (2018) J. Fischer _et al._ , Phys. Rev. B 97, 014417 (2018).
  * Lebrun  _et al._ (2019) R. Lebrun _et al._ , Communications Physics 2, 50 (2019).
  * Jia  _et al._ (2011) X. Jia, K. Liu, K. Xia, and G. E. W. Bauer, Europhys. Lett. 96, 17005 (2011).
  * Chen  _et al._ (2013) Y.-T. Chen _et al._ , Phys. Rev. B 87, 144411 (2013).
  * Vlietstra  _et al._ (2014) N. Vlietstra _et al._ , Phys. Rev. B 90, 174436 (2014).
  * Schwarze  _et al._ (2015) T. Schwarze _et al._ , Nat. Mat. 14, 478–483 (2015).
  * Hanke  _et al._ (2020) J.-P. Hanke _et al._ , Phys. Rev. B 101, 014428 (2020).
  * Keesman  _et al._ (2015) R. Keesman _et al._ , Phys. Rev. B 92, 134405 (2015).