# SNR Spectra as a Quantitative Model for Image Quality in Polychromatic X-Ray
Imaging

M. Ullherr  Fraunhofer Development Center X-ray Technology EZRT, Josef-Martin
Weg 63, 97074 Würzburg, Germany  Corresponding author:
<EMAIL_ADDRESS>S. Zabler  Fraunhofer Development
Center X-ray Technology EZRT, Josef-Martin Weg 63, 97074 Würzburg, Germany

###### Abstract

In polychromatic x-ray imaging for nondestructive testing, material science or
medical applications, image quality is usually a problem of detecting sample
structure in noisy data. This problem is typically stated this way: As many
photons as possible need to be detected to get a good image quality.

We instead propose to use the concept of signal detection, which is more
universal. In signal detection, it is the sample properties which are
detected. Photons play the role of information carriers for the signal. Signal
detection for example allows modeling the effects which polychromaticity has
on image quality.

$\mathit{SNR}$ spectra (= spatial $\mathit{SNR}$) are used as a quantity to
describe if reliable signal detection is possible. They include modulation
transfer and phase contrast in addition to noisiness effects. $\mathit{SNR}$
spectra can also be directly measured, which means that theoretical
predictions can easily be tested.

We investigate the effects of signal and noise superposition on the
$\mathit{SNR}$ spectrum and show how selectively not detecting photons can
increase the image quality.

##  1 Introduction

X-ray imaging light sources are relatively weak light sources, so a detection
process in x-ray imaging is often noise-limited. This is especially true for
computed tomography (CT), due to the need to acquire many images and the
correspondingly longer measurement times. There are other detection problems
in x-ray imaging that are not noise-limited, which we do not consider here.

It is therefore a natural assumption that, of the photons which pass through
the sample, the detected fraction needs to be as high as possible to achieve a
high image quality. Also, there should be as many photons generated as
possible. This can lead to an understanding of image quality where it is
implicitly assumed that detecting photons is the main aspect. As this sets the
focus on detecting photons, possibly at the exclusion of other considerations,
we will call this concept ”photon detection”. Examples for work which
(implicitly) use this concept to evaluate image quality are [1, 2, 3, 4, 5]

While photon detection is a sufficient model in some cases, it is an
oversimplification in others: For example if polychromaticity or modulation
transfer are important effects in x-ray imaging. Here polychromaticity means
that different x-ray energies contribute to a single image. Generally, if
different photons in one image have very different properties, photon
detection is an oversimplification. To fully model these cases, we
fundamentally consider the problem of signal detection in imaging, starting
with a problem definition:

Detection is a measurement process in which a specific signal is reliably
differentiated from noise.
---

To make this definition complete, the terms used are defined and explained in
table 1.

term |  definition
---|---
signal |  Sample properties and structure which are intended to be measured.
data |  Result of a measurement.
noise |  Deviation between measured values and signal.
image quality |  Degree of reliability with which the sample structure of interest can be detected in a measured image.
modulation transfer |  Image blur (signal deterioration) caused by the detection apparatus (source size, diffraction limit, …). Quantified by the modulation transfer function (MTF).
poly-chromaticity |  Photons have a broad energy distribution.
Table 1: Definitions for basic terms in the signal detection model. While it
is necessary and useful to define the truth (signal) as a theoretical concept,
it is impossible to know exactly. Using photons as information carriers for
the signal means that there is always (at least) Poisson noise.

The main difference to the concept of photon detection is that we consider the
detection of a signal, not the detection of photons. Therefore we call this
concept ”signal detection”. Photons are used to probe the signal, they are
information carriers, not the signal itself.

In the case of x-ray imaging, the signal is the spatial distribution
(structure) of the x-ray photon interaction strengths of the sample. The
interaction can either be absorption or phase effects, which respectively
correspond to the imaginary and real parts of the index of refraction.
Usually, only parts of the sample structure is of interest. We define the term
”image quality” to describe how well these can be detected.

This work is split into two parts. In the first part (section 2), we will take
a closer look at what a quantitative model for image quality requires to be
sufficiently accurate. We then consider if the common approaches satisfy these
requirements and explain why SNR spectra are a good quantitative model. In the
second part (the rest of this work), we develop a framework for understanding
and evaluating polychromatic image quality based on SNR spectra.

##  2 Image Quality Measures

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 1: Simulated projection images of randomly placed balls of different
sizes with the same number of balls for each size (diameters: 6.0, 4.2, 3.0,
2.1, 1.5). Images b) and c) are generated by adding another image to d). For
g)–i), $\sigma$ is the standard deviation of a Gaussian MTF and 1000 counts
were used.

###  2.1 Effects to Include

We now want to find a physical model to describe image quality in signal
detection quantitatively. One of the basic difficulties when designing a
physical model is deciding which effects to include and which effects to
ignore e.g. by an approximation. A careful consideration of appropriate
approximations is required, otherwise the physical model will not describe
reality correctly in the relevant cases.

The physical model considered in this work is that of a numerical value to
describe image quality in a signal detection context. As a value to describe a
quantity, we will call these different values ”measures”.111 Note that these
types of measures do not satisfy additivity – as would be expected from the
use in mathematics. These measures are then used to optimize an imaging device
or answering the question if one measurement technique is better than another.
If the measure does not correctly describe the quantity, the wrong case might
be considered superior.

Figure 2: Energy dependency of the signal strength for x-ray absorption for
some materials. A variation by two orders of magnitude is realistic for broad
polychromatic spectra. Approximating $\mu$ as constant is generally only
appropriate for monochromatic or very high energy imaging (see Al curve for
$>$ 60 keV or Fe for $>$ 150 keV). This is caused by the fact that Compton
scattering becomes the stronger interaction at higher energies and its energy-
dependency is weak.

In polychromatic x-ray imaging, effects from the following categories need to
be included:

  * (A)

Noisiness (intensity/Poisson noise and other)

  * (B)

Modulation transfer (MTF)

  * (C)

Signal strength (absorption, phase or other)

  * (D)

Superposition of contributions with different A—C

  * (E)

Presence of artifacts ($\approx$ noise)

Example images for how these effects influence image quality are given in Fig.
1. They consist of randomly placed balls of different size but identical
maximal absorption length. These balls are then blurred with a MTF and Poisson
noise is simulated. The examples b) and c) are both superpositions of d) and
another image with low image quality. In both cases, the resulting image has a
lower image quality than d) alone. It can be seen that all effects have the
same result: Small objects can or cannot be reliably detected depending on
different properties of the imaging setup. For example noisiness (A) and MTF
(B) have the same resulting effect and should therefore be quantified by a
unified measure.

Note that from Fig. 1 d)–f) we can see that small differences in the image
quality are not noticeable to a human observer. A difference by a factor of
three can be differentiated well, much smaller differences would not be. This
makes image quality optimizations by human estimate very imprecise.

The effects A-C scale the image quality. Ignoring one of these effects can
give rise to severe misinterpretations. Effects D are important if one of the
effects A, B or C differs greatly for different detected photons, even if the
effect themselves would not need to be considered. In Fig. 2, an example for a
strong energy dependency of the signal strength is shown. Different
intensities (A) and varying signal strengths (C) have the same effect in a
single image, but different effects in a superposition. Otherwise, both
effects could be modeled as one effect. Effects E are included for
completeness, but are not considered further in this work.

The two commonly used measures to describe and optimize image quality are
$\mathit{DQE}$ (= photon detection efficiency) and $\mathit{SNR}$ (=
ability/probability to detect signal). Both are used in temporal (scalar) and
in spatial (frequency-dependent) forms.

###  2.2 Signal to Noise Ratio

The signal to noise ratio gives the strength of the signal in relation to the
strength of the noise. It can be understood as an (inverse) relative
measurement error. This concept is based on the fact that signal can only be
reliably detected if it is stronger than the noise (relative error $\ll 1$).
While it is clear that higher $\mathit{SNR}$ is better, defining a
sufficiently high $\mathit{SNR}$ depends on many factors and can not be done
generally. The definitions of the temporal $\mathit{SNR}_{t}$ and the
$\mathit{SNR}$ spectrum are:

$\displaystyle\mathit{SNR}_{t}$ $\displaystyle=\frac{\text{pixel\ mean\
signal}}{\text{pixel\ standard\ deviation}}\big{(}=\mathit{CNR}_{t}\big{)}$
(1) $\displaystyle\mathit{SNR}(u)$ $\displaystyle=\frac{\text{image\ signal\
power\ spectrum}}{\text{image\ noise\ power\ spectrum}}$ (2)

$u$ is the spatial frequency and $\mathit{SNR}(u)$ is the $\mathit{SNR}$
spectrum. ”Spectrum” here means image spatial frequency spectrum, not light
wavelength spectrum. Note that by ”signal”, we always mean sample structure,
not the light intensity. Our definition of $\mathit{SNR}_{t}$ is identical to
a contrast to noise ratio ($\mathit{CNR}_{t}$) – the contrast _is_ the signal.
This is consistent with how $\mathit{SNR}(u)$ is defined.

While $\mathit{SNR}_{t}$ is suitable to describe the lower noisiness of images
with longer integration times, comparing cases with differences in the
measurement setup (MTF differs e.g. due to a different/hardened x-ray
spectrum, different screen thickness, …) potentially leads to false
comparisons. Optimizing $\mathit{SNR}_{t}$ may thus lead to a lower image
quality if this comes at the cost of a worse MTF.

Determining and evaluating both $\mathit{SNR}_{t}$ and MTF can show this
problem. It does not answer the question as to which case is optimal, which
the $\mathit{SNR}(u)$ does.

$\mathit{SNR}_{t}$ is a linear quantity while $\mathit{SNR}(u)$ is a squared
quantity. The former is proportional to the signal amplitude while the latter
is proportional to measurement integration time.

###  2.3 Detective Quantum Efficiency

A $\mathit{DQE}$ is defined as a $\mathit{SNR}$ transfer function to describe
and optimize one part of the imaging device. Usually it is applied to describe
the x-ray detector. It is generally calculated in the form

$\displaystyle\mathit{DQE}=\frac{\mathit{SNR}_{\text{out}}}{\mathit{SNR}_{\text{in}}}=\frac{\mathit{SNR}_{\text{detected}}}{\mathit{SNR}_{\text{ideal}}}\in[0,1[$
(3)

This is intended to give an absolute optimum if $\mathit{DQE}=1$. There exist
many derived expressions for certain applications (simulation [2, 6] or
measurement[7, 8]). The derivations assume monochromaticity but are applied to
polychromatic imaging anyway. They assume that signal strengths cancel out and
detector properties can be averaged over the x-ray spectrum. The common
simplified expression for the $\mathit{DQE}$ of an indirect detector is:

$\mathit{DQE}(u)=\frac{a}{1+c^{-1}H_{v}(u)^{-2}}$ (4)

where $a$, $c$ and $H_{v}(u)$ are the spectrally averaged x-ray absorption
efficiency, x-ray photon conversion factor and optical (scattering) MTF.222If
$a$ is included varies, but it must be included when considering the whole
x-ray photon detection process. Effects from different samples are excluded by
design in such a simplification, which means that effects (C) and (D) are
neglected.

$\mathit{SNR}_{\text{ideal}}$ and $\mathit{SNR}_{\text{detected}}$ depend on
the sample, which has a x-ray energy dependent transmission. This effective
x-ray spectrum would need to be considered to correctly compute a
polychromatic $\mathit{DQE}$ and can be strongly sample-dependent.

Using a monochromatic $\mathit{DQE}$ at different energies does not pose such
problems. This is the standard approach for characterizing detectors for
optical light (CMOS/CCD) [9], in which a similar quantity is called quantum
efficiency (QE). It is the only application of the $\mathit{DQE}$ concept in a
polychromatic context which is accurate.

###  2.4 Indirect Detection

Generally, two types of x-ray imaging detectors can be differentiated:

  * $\circ$

Direct detectors: These detectors count single x-ray photon detection events.
They are principally able to determine the x-ray energy of a detected photon.
They are also called counting detectors.

  * $\circ$

Indirect detectors: These detectors first convert x-ray photons into visible
light which is then detected. They usually integrate all events in a certain
time frame, weighted by the energy of the x-ray photon. They are also called
(energy) integrating detectors.

Due to the fact that indirect detectors are currently much cheaper to produce,
most x-ray imaging setups use this kind of detector.

For indirect detection, the noise power spectrum $N(u)$ is not a white
spectrum, but is instead influenced by the detector MTF in the following way
[10]:

$\displaystyle N(u,E)=I(E)\left[c(E)^{2}H_{v}(u,E)^{2}+c(E)\right]$ (5)

where $E$ is the x-ray energy, $H_{v}(u,E)$ is the optical detector MTF,
$I(E)$ is the x-ray intensity and $c(E)$ is the conversion factor for x-rays
to visible light. This equation is only valid for monochromatic x-rays and for
a simplified imaging setup, but the general shape of $N(u)$ is identical for
polychromatic spectra on real setups. In a CT measurement, the reconstruction
algorithm will also influence volume image noise.

The pixel noise $\sigma_{t}$ (temporal standard deviation) for a noise power
spectrum $N(u)$ of a projection is given by Parseval’s relation [11]:

$\sigma_{t}^{2}=\frac{1}{A}\sum_{x,y}N(u_{x,y})$ (6)

where $N(u_{x,y})$ is one discrete Fourier coefficient and $A$ is the number
of these coefficients. Computing or measuring a $\sigma_{t}$ for an indirect
detector does not give useful results even for noisiness effects, because
noise at all frequencies is averaged, although it may not decrease image
quality at the structure size that is of interest.

If the image quality is reduced by decreasing $H_{v}(u,E)$, both
$\mathit{SNR}_{t}$ and $\mathit{DQE}_{t}$ will increase due to a decreasing
$\sigma_{t}$. Changing $H_{v}(u,E)=1$ to $H_{v}(u,E)=0$ decreases $\sigma_{t}$
by a factor of $\sqrt{c+1}$. If technically feasible, $c>100$ is achieved in
detectors built today. This effect can therefore be large. In other words,
sharper indirect detectors are also temporally noisier. The $\mathit{SNR}$
spectrum instead indicates that they have relatively weaker noise, which is
correct.

Similar effects are produced by artifacts. Beam hardening artifacts or ring
artifacts may e.g. decrease the $\mathit{SNR}(u)$ at lower frequencies but
fine detail can still be detected well. A temporal $\mathit{SNR}$ cannot
differentiate between low-frequency and high-frequency noise.

###  2.5 SNR Spectra

| effect included
---|---
| A | B | C | D | E
$\mathit{DQE}_{t}$ | yes∗ | no | no | no | no
$\mathit{DQE}(u)$ | yes∗ | yes∗ | no | no | no
$\mathit{SNR}_{t}$ | yes∗ | no | yes | yes∗ | no
$\mathit{SNR}(u)$ | yes | yes | yes | yes | no
Table 2: Summary of image quality measures and represented physical effects.
Entries marked with a $*$ include assumptions or approximations that are
incorrect in common cases.

Which of the image quality measures discussed takes which physical effect into
account is summarized in table 2.

Of those measures considered, only $\mathit{SNR}$ spectra are left as a
candidate for a generally useful image quality measure – they do not suffer
from any of the problems discussed previously. While the other image quality
measures do work in special cases, their restrictions have often been ignored
and can lead to severe misinterpretations of physical effects on image
quality.

##  3 Signal Detection Model

Understanding image quality in the context of signal detection requires that
we model polychromatic effects. To do so, we will use $\mathit{SNR}$ spectra.
These are computed from the fraction of the signal power spectrum to the noise
power spectrum, see eq. (2). Signal detection image quality can therefore be
reduced to the properties of these power spectra, which we will consider in
the following.

###  3.1 Classes of Photons

To understand how the different detection properties (e.g. due to
polychromaticity) affect $\mathit{SNR}$ spectra, we will first consider the
signal and noise spectra of photons with the same detection properties. This
may be e.g. a monochromatic thin-screen model which is then used as a building
block for the polychromatic thick-screen model.

Following the use in set theory, we will call this a class of photons: The
(detected) x-ray photons that have (approximately) identical detection
properties. The two most important variables for detection properties are $E$,
the x-ray energy and $z$, the position within the scintillation screen. The
MTF depends on the latter and the depth intensity distribution is energy-
dependent.

The actual image then consists of contributions from classes of photons – the
superposition of the images of all individual classes.

###  3.2 Signal and Noise Superposition

To model $\mathit{SNR}$ spectra of polychromatic images, the power spectra of
signal and noise must be calculated for a superposition of images. In the
strictest sense, this can be seen as a superposition of x-ray detection events
(one image = one event). We thus need to derive the superposition equations
for signal and noise power spectra.

If $d_{\Sigma}(x)$ is the polychromatic superimposed intensity image at the
pixel position $x$ and $d(x,E)$ the image contribution at energy $E$ (= photon
class), we get:

$d_{\Sigma}(x)=\sum_{E}d(x,E)$ (7)

where the $\sum$ may also represent an energy integral, in which case $d(x,E)$
is a density. In general, this sum may also be a multiple sum/integral over
several variables ($E$, $z$, …), where one combination of variables
corresponds to a photon class.

The linear additive noise model can be used to split $d$ into signal $s$ and
noise $n$:

$\displaystyle d(x)$ $\displaystyle=s(x)+n(x)$ (8)

The fact that there is a transformation between physical signal (volume
distribution) and $s(x)$ (projection image) is omitted here because it does
not influence the results. In a volume $\mathit{SNR}$ analysis, the noise is
also transformed by the CT reconstruction. We can derive polychromatic
expressions for signal and noise by adding up intensities in image space:

$\displaystyle s_{\Sigma}(x)$ $\displaystyle=\sum_{E}s(x,E)$ (9)
$\displaystyle n_{\Sigma}(x)$ $\displaystyle=\sum_{E}n(x,E)$ (10)

The image power spectrum is calculated as:333In the following we denote power
spectra by capital letters, not the Fourier transformed quantities (except for
the MTF $H$). Note that if power spectra are measured, corresponding
measurement errors occur.

$\displaystyle D_{\Sigma}(u)$
$\displaystyle=\left|\mathcal{F}\left\\{d_{\Sigma}(x)\right\\}(u)\right|^{2}=S_{\Sigma}(u)+N_{\Sigma}(u)$
(11)

because noise is uncorrelated with any other signal/noise and the
corresponding mixed terms vanish. Similarly, for the noise power spectrum
itself we get:

$\displaystyle N_{\Sigma}(u)$
$\displaystyle=N_{0}(u)+\left|\mathcal{F}\left\\{\sum_{E}n(x,E)\right\\}(u)\right|^{2}$
(12) $\displaystyle=N_{0}(u)+\sum_{E}N(u,E)$ (13)

where $N_{0}$ are noise contributions in addition to photon noise, e.g. camera
readout/dark noise. Scattered photons are counted in the noise sum (energy of
the scattered photon), but usually do not contribute to the signal.

For the signal, we get:

$\displaystyle S_{\Sigma}(u)$
$\displaystyle=\left|\mathcal{F}\left\\{\sum_{E}s(x,E)\right\\}(u)\right|^{2},$
(14)

which is simplified in the next section.

###  3.3 Local Area Approximation

If we consider the image quality in a small area (locally), we can approximate
the intensity as constant within that area. In the following, intensity thus
stands for the number of x-ray photons (units of 1) detected in a local area.
A local area is much larger than a single pixel and contains weak sample
structure. This also allows us to assume that for absorption, the signal
strength is proportional to the absorption strength. The proportionality is
given by a linear approximation of the absorption curve (Lambert-Beer law).

Because in the end we are interested in finding an optimum of the unit-free
quantity $\mathit{SNR}$, physical units are not considered. They can be
restored if needed. We can thus write a simplified model for the signal and
its power spectrum as:

$\displaystyle s(x,E)$
$\displaystyle=I(E)\left[\kappa(E)\rho(x)\right]\,*\,h(x,E)$ (15)
$\displaystyle S(u,E)$ $\displaystyle=P(u)\
\big{[}I(E)\kappa(E)H(u,E)\big{]}^{2}$ (16)

The signal is a product of the detected intensity $I(E)$, an energy
independent spatial matter density $\rho(x)\leftrightarrow P(u)$, the MTF
$H(u,E)$ and an interaction strength $\kappa(x,E)$. We have used
$\kappa(x,E)=\kappa(E)$ (single material sample) for simplicity—considering
multi material samples would not give new effects. Note that we assume the
imaging device to be a linear translationally invariant (LTI) system in the
local area.

$\kappa(E)$ can model almost any signal generating physical effect. This is
usually x-ray absorption or phase contrast. In an cases like imaging with
x-ray optics or grating interferometry [12], $\kappa(x,E)$ is also determined
by the energy band for which the x-ray optics or gratings are optimized.

For a signal sum eq. (14) becomes:

$\displaystyle S_{\Sigma}(u)$
$\displaystyle=\mathrm{P}(u)\left[\sum_{E}I(E)\kappa(E)H(u,E)\right]^{2}$ (17)

The influence of the sample structure $P(u)$ is a global scaling factor. The
right hand side originally contains the square of the absolute value, but all
sum contributions are real valued. It can be easily seen that signal sums are
difficult to evaluate, because many physical effects need to be
considered—even in a basic model.

Signal power spectra of large structures show higher amplitude at lower
spatial frequencies and $P(u)$ is then usually monotonically decreasing
towards higher frequencies. The reason is that sample structures are long-
range correlations in the signal, and the Wiener-Chinchin-Theorem translates
this to the power spectrum shape. Larger structural details therefore have an
intrinsically higher $\mathit{SNR}$ spectrum value and are easier to detect. A
lower $\mathit{SNR}(u)$ usually has the effect that the detail resolution of a
measurement gets worse.

It is possible to define a quantity that is not dependent on the amount of
structure in the sample or on the shape of the signal power spectrum. We will
call this the detection effectiveness spectrum, which is derived from eq. (13)
and eq. (17) and defined as:

$\displaystyle\mathit{DE}(u)$
$\displaystyle=\frac{\mathit{SNR}_{\Sigma}(u)}{P(u)\eta}=\frac{S_{\Sigma}(u)}{N_{\Sigma}(u)P(u)\eta}$
$\displaystyle=\frac{\big{[}\sum_{E}I(E)\kappa(E)H(u,E)\big{]}^{2}}{\big{[}N_{0}(u)+\sum_{E}N(u,E)\big{]}\eta}$
(18)

Introducing the cost function $\eta$ further allows to compare measurements
with e.g. different measurement times ($\eta$ = time). In medical imaging,
using the dose as a cost function can be appropriate. Optimizing the
$\mathit{SNR}(u)$ will give the same result as optimizing the $\mathit{DE}(u)$
if $\eta$ is constant. In the following considerations, the $\mathit{SNR}(u)$
could usually be replaced by the $\mathit{DE}(u)$.

The $\mathit{DE}(u)$ does not depend on the specific sample structure and thus
describes only the imaging setup itself. Its main advantage compared to the
$\mathit{SNR}(u)$ is that the $\mathit{DE}(u)$ allows comparable measurements.
In contrast to e.g. a $\mathit{DQE}$, this quantity models polychromatic
effects. Sample-dependent properties of an imaging setup are included by the
$\mathit{DE}(u)$.

###  3.4 SNR Spectra of Superpositions

One important aspect of $\mathit{SNR}$ spectra are the effects that come from
superimposing images with different properties, e.g. different $\kappa$ or
MTF. From eq. (13) and eq. (17) we get:

$\mathit{SNR}_{\Sigma}(u)=P(u)\frac{\big{[}\sum_{k}I_{k}\kappa_{k}H_{k}(u)\big{]}^{2}}{N_{0}(u)+\sum_{k}N_{k}(u)}$
(19)

where $k$ is the index for different contributions to the resulting image – a
contribution may consist of a single photon class or can itself be a sum of
photon classes. While both sums are linear superpositions, the square in the
numerator makes many simplifications impossible. An important property of eq.
(19) is that the signal contributions are weighted with the signal strength
$\kappa$ and MTF $H$, while the noise contributions are not.

We can see from eq. (19) that adding two images with greatly differing
$\mathit{SNR}$ (due to relative differences in signal strength or MTF) will
result in an image with a lower $\mathit{SNR}$ than that of the one image with
the higher $\mathit{SNR}$. The following triangle inequality holds:

$\mathit{SNR}_{1+2}(u)\leq\mathit{SNR}_{1}(u)+\mathit{SNR}_{2}(u)$ (20)

Detecting more photons (image 2) in addition to the reference situation (image
1) increases the overall $\mathit{SNR}$ if:

$\displaystyle\mathit{SNR}_{1+2}(u)$ $\displaystyle>\mathit{SNR}_{1}(u)$ (21)
$\displaystyle\Leftrightarrow\frac{\big{[}I_{1}\kappa_{1}H_{1}+I_{2}\kappa_{2}H_{2}\big{]}^{2}}{N_{0}+N_{1}+N_{2}}$
$\displaystyle>\frac{\big{[}I_{1}\kappa_{1}H_{1}\big{]}^{2}}{N_{0}+N_{1}}$
(22)
$\displaystyle\Leftrightarrow\frac{I_{1}\kappa_{1}H_{1}+I_{2}\kappa_{2}H_{2}}{I_{1}\kappa_{1}H_{1}}$
$\displaystyle>\sqrt{\frac{N_{0}+N_{1}+N_{2}}{N_{0}+N_{1}}}$ (23)
$\displaystyle\Leftrightarrow\frac{I_{2}\kappa_{2}H_{2}}{I_{1}\kappa_{1}H_{1}}$
$\displaystyle>\sqrt{1+\frac{N_{2}}{N_{0}+N_{1}}}-1$ (24)

If this condition is violated, the additional photons (image 2) would decrease
$\mathit{SNR}$ if they are detected. In this equation, the absolute signal
strengths are not important, but the relative strengths are.

For an infinitely small addition ($N_{2}\ll N_{1}$) we can approximate the
square root above and get:

$\displaystyle\frac{I_{2}\kappa_{2}H_{2}}{N_{2}}$
$\displaystyle>\frac{1}{2}\frac{I_{1}\kappa_{1}H_{1}}{N_{1}+N_{0}}$ (25)

While this is an interesting limit case, it does not have any practical
application. We define these naturally arising quantities as the ”signal
detection strength”:

$\displaystyle\mathit{SDS}(u,E)$
$\displaystyle=\frac{I(E)\kappa(E)H(u,E)}{N(u,E)}$ (26)

A $\mathit{SDS}$ can be calculated either for a photon class or for a set of
photon classes ${\\{E\\}}$:

$\mathit{SDS}(u,\\{E\\})=\frac{\sum_{\\{E\\}}I(E)\kappa(E)H(u,E)}{\sum_{\\{E\\}}N(u,E)}$
(27)

For $N_{0}\ll N_{1}$ we can rewrite eq. (25) to:

$\mathit{SDS}_{2}(u,E)>\frac{1}{2}\mathit{SDS}_{1}(u,E)$ (28)

For an ideal detector ($I=N$), the $\mathit{SDS}$ simplifies to:

$\mathit{SDS}(u,E)=\kappa(E)H(u,E)$ (29)

Note that for absorption imaging, $H(u,E)\leq 1$ and
$\mathit{SDS}(u,E)\leq\kappa(E)$. For phase contrast imaging, $H(u,E)$ may be
$>1$ (or even $\gg 1$).

From eq. (28) we can see that a small intensity addition requires at least
half the $\mathit{SDS}$ of the existing image to contribute positively. On the
other hand, for a large intensity addition ($I_{2}\gg I_{1}$),
$\mathit{SNR}_{1+2}>\mathit{SNR}_{1}$ is always true. As an example for an
intermediate case, if $N$, $H$ and $I$ are equal for both images, the
condition becomes $\kappa_{2}>0.41\kappa_{1}$.

###  3.5 Weighted Superposition

As discussed in [13], applying an energy-dependent weighting factor before the
superposition of the image contributions can increase the $\mathit{SNR}$. This
is called ”energy weighting” and there are two variants:

  * $\circ$

Multiplication with an energy-dependent weighting factor after the detection
of the x-ray photons. This is what is usually discussed as energy weighting.
Because sufficiently good energy-resolving 2D imaging detectors are currently
not technically feasible, this is mostly a theoretical concept.
We will call this effect ”computational energy weighting” (CEW).
A CEW weight $w$ scales $\mathit{SDS}$ with $1/w$.

  * $\circ$

Implicit weighting by the different detected intensities $I(E)$ (effective
x-ray spectrum). This is influenced e.g. by the source spectrum or the
detector absorption efficiency spectrum. This weighting can be changed by
purposefully decreasing the detected intensity e.g. in specific energy ranges.
We will call this effect ”detection energy weighting” (DEW).

In both cases, optimizing the weighting can lead to an increase in
$\mathit{SNR}$. The increases vary for different experimental conditions
between a few percent and factors $>10$. The larger effect sizes are for very
broad spectra that are not useful without energy weighting.

For the CEW, if we multiply with a weight $w(u,E)$, we get:

$\mathit{SNR}_{\text{CEW}}(u,w)=\mathrm{P(u)}\frac{\big{[}\sum
w(u,E)I(E)\kappa(E)H(u,E)\big{]}^{2}}{\sum w(u,E)^{2}N(u,E)}$ (30)

To compute the optimal energy weighting function $w_{\text{opt.}}$ for CEW, we
determine the maximum of the $\mathit{SNR}$ spectrum depending on $w(u,E)$:

$\frac{\mathrm{d}\mathit{SNR}_{\text{CEW}}(u,w)}{\mathrm{d}w}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{!}}}}{{=}}}0$
(31)

which is solved by any CEW weight proportional to

$w_{\text{opt.}}(u,E)\propto\frac{I(E)\kappa(E)H(u,E)}{N(u,E)}=\mathit{SDS}(u,E)$
(32)

This is an extension of the known result
$w_{\text{opt.}}(u,E)\propto\kappa(E)$ [13]. For a detector with a limited
energy resolution, eq. (27) can be used. To avoid artificial image blurring,
it is sufficient to only apply a relative blurring factor in this way:

$w_{\text{opt.}}(u,E)=\frac{I(E)\kappa(E)H(u,E)}{N(u,E)H_{0}(u)}$ (33)

where $H_{0}(u)$ can e.g. be the energy-averaged MTF. Note that the
$\mathit{SNR}$ spectrum is not changed when a Fourier filter is applied, but
the image may be easier or more difficult to interpret. After application of
the optimal CEW, the signal detection strength $\mathit{SDS}$ is constant with
respect to $E$:

$\mathit{SDS}_{\text{opt.\ CEW}}(u,E)=\text{const}(E)$ (34)

The special case of a task-independent (t.i.) measurement case is defined as a
case with energy independent signal detection strength (without weighting):

$\displaystyle\mathit{SDS}_{\text{t.i.}}(u,E)$ $\displaystyle=\text{const}(E)$
(35) $\displaystyle\Rightarrow w_{\text{opt}|\text{t.i.}}(u,E)$
$\displaystyle=\text{const}(E)$ (36)
$\displaystyle\Rightarrow\mathit{SNR}_{\Sigma|\text{t.i.}}(u)$
$\displaystyle=\sum_{E}\mathit{SNR}_{\text{t.i.}}(u,E)$ (37)

In this case, the optimal $\mathit{DQE}$ also gives the optimal
$\mathit{SNR}$. This is the simplest possible case which is often implicitly
used to understand detection efficiency. To do so, eq. (35) must be assumed to
hold true approximately. This is only appropriate in some special cases, e.g.
monochromatic imaging.

If we compare eq. (34) and eq. (35), we can see that applying an optimal CEW
has the effect that the result of the CEW appears to be task-independent. The
CEW itself is of course task-dependent. This also means that task-independency
is the case in which there is an optimal superposition of the $\mathit{SNR}$
contributions.

For the optimal DEW, there appears to be no analytical expression. The optimal
DEW is thus the case with the maximal $\mathit{SNR}$. Eq. (24) can be used to
find approaches that may increase $\mathit{SNR}$. The general principle is
that if the x-ray photons at specific energy ranges do not contribute
positively to $\mathit{SNR}$, a reduction of their intensity increases
$\mathit{SNR}$.

##  4 Examples

###  4.1 Two Images Superposition

description | $\kappa_{2}$ | $\mathit{SNR}_{2}$ | $\mathit{SNR}_{1+2}$
---|---|---|---
a) photoabsorption | 0.14 | 0.020 | 0.65
b) phase contrast | 0.25 | 0.063 | 0.78
c) no gain, eq. (24) | 0.41 | 0.17 | 1
d) task-independent | 1 | 1 | 2
e) optimal CEW for a) | 0.14 | 1 | 1.02
f) optimal DEW for a) | 0.14 | 1 | 1
Table 3: Examples for the $\mathit{SNR}(u_{a})$ of image sums for different
physical effects represented by $\kappa$ with $\kappa_{1}$ = 1 and
$\mathit{SNR}_{1}$ = 1. Cases a) and b) are e.g. iron for 30 keV and 60 keV.
Note that for a) to d), the fraction of detected photons ($\propto$
polychromatic $\mathit{DQE}$) is twice as high for the sum image, but for a)
and b) the image quality is lower. See section 3.5 for cases e)$+$f).

Figure 3: Examples for the dependency of the $\mathit{SNR}$ on the weighting
factor for CEW (top) and DEW (bottom). The weighting factor factor of image 2
is $w_{2}$. The annotations correspond to table 3. For DEW, the maximal
$\mathit{SNR}$ is either at $I_{2}=0$ or at finite $I_{2}$, depending on the
maximal value of $I_{2}$. The edge case $\kappa_{2}=0.5\kappa_{1}$ is the
lowest possible $\kappa_{2}$ where $\mathit{SNR}_{1+2}\geq\mathit{SNR}_{1}$
for all $I_{2}$.

As the simplest possible example for $\mathit{SNR}$ spectra of an image
superposition, we will consider adding two monochromatic images:

$d_{1+2}(E)=d_{1}\delta(E-E_{1})+d_{2}\delta(E-E_{2})$ (38)

We assume the following conditions:

$\displaystyle P(u_{a})H_{1}(u_{a})^{2}=P(u_{a})H_{2}(u_{a})^{2}=1;$ (39)
$\displaystyle I_{1}=I_{2}=1;\ \ E_{1}/E_{2}=0.5$ (40)

for a structure size of $u_{a}$. This effectively excludes modulation transfer
effects from our examples—this can be done in a theoretic model without loss
of generality,444 Including $H(u_{a})$ as a factor in $\kappa$ restores the
full properties. but would be almost useless for a real application.

This is a realistically polychromatic case whose effect sizes approximate real
cases. Inserting all these values into equations (13) and (17) yields:

$\displaystyle N_{1+2}(u_{a})$ $\displaystyle=I_{1}+I_{2}$ (41) $\displaystyle
S_{1+2}(u_{a})$ $\displaystyle=(I_{1}\kappa_{1}+I_{2}\kappa_{2})^{2}$ (42)
$\displaystyle N_{1}(u_{a})$ $\displaystyle=I_{1};$ $\displaystyle
N_{2}(u_{a})$ $\displaystyle=I_{2}$ (43) $\displaystyle S_{1}(u_{a})$
$\displaystyle=(I_{1}\kappa_{1})^{2};$ $\displaystyle S_{2}(u_{a})$
$\displaystyle=(I_{2}\kappa_{2})^{2}$ (44)

where we assume an ideal photon counting detector. Using the specific values
listed above gives:

$\displaystyle\mathit{SNR}_{1}$ $\displaystyle=\mathit{SNR}_{1}(u_{a})=1$ (45)
$\displaystyle\mathit{SNR}_{2}$ $\displaystyle=\kappa_{2}^{\ 2}$ (46)
$\displaystyle\mathit{SNR}_{1+2}$ $\displaystyle=\frac{(1+\kappa_{2})^{2}}{2}$
(47)

For some examples of $\kappa_{2}$ the $\mathit{SNR}$ values for image sums are
computed in table 3. From these examples we can see that the additionally
detected intensity $I_{2}$ may reduce the overall image quality significantly
if the $\mathit{SNR}_{2}$ of this additional intensity is low compared to the
$\mathit{SNR}_{1}$ of $I_{1}$. A difference in a MTF may also produce this
effect.

The reason for this behavior is easily explained by the fact that if
$\kappa_{1}^{2}\gg\kappa_{2}^{2}$, $I_{2}$ contributes little signal but all
of its noise to the image sum. A lower detected intensity at $E_{2}$ would
decrease the noise.

For this example, a high enough intensity can compensate for a low signal
strength in the following way:

$\mathit{SNR}_{1+2}\geq\mathit{SNR}_{1}\Leftrightarrow I_{2}\geq
I_{1}\frac{1-2\kappa_{2}}{\kappa_{2}^{\ 2}}$ (48)

This simple example allows us to evaluate the effect (non-optimal) energy
weighting has on the $\mathit{SNR}$, which is shown in Fig. 3 for CEW and DEW.
The examples only qualitatively represent realistic cases, but demonstrate the
effects which different energy weightings have on $\mathit{SNR}$.

We can see that for CEW, there is always a unique maximum at
$w_{2}=\kappa_{2}$. For DEW, there is a unique minimum, while the maximum is
obtained for $I_{2}\rightarrow\infty$, and for $\kappa_{2}<\kappa_{1}/2$ an
additional local maximum exists at $I_{2}=0$. If we assume that a maximal
value for $I_{2}$ is given by the physical circumstances (e.g. source
spectrum), the optimal weighting is either to use this maximal value or use
$I_{2}=0$ (depending on $\kappa_{2}$). The optimal DEW weight can thus be
interpreted as a mask function that is either 0 or 1. This is usually a step
function that is 1 at low energies and 0 above some threshold.

###  4.2 Simple Polychromatic Simulations

Figure 4: Simulated (photo-)absorption $\mathit{SNR}(u_{a})$ (top) and x-ray
inline phase contrast $\mathit{SNR}(u_{a})$ (bottom) curves for a
bremsstrahlung source with constant source power. The corresponding detected
x-ray spectra are shown in the middle, where the highest energy in a curve is
the tube voltage.

$\mathit{SNR}$ simulations for simple x-ray spectra are given in Fig. 4. They
are generated by analytically computing $\mathit{SNR}(u)$, using eq. (19).
Both examples use the same source spectra, detector and sample, but different
energy-dependencies of the signal strength. These examples demonstrate the
magnitude of the effect different energy weighting effects have on
$\mathit{SNR}$, they are not meant to describe an actual setup. See Fig. 2 for
the absorption coefficient for aluminium (Al), which gives the energy-
dependent sample transparency ($0.1\text{\,}\mathrm{mm}$ thickness) and the
signal strength for the absorption signal (approximately $\kappa(E)\propto
E^{-3}$). The x-ray inline phase contrast signal strength is assumed to be
$\kappa(E)\propto E^{-2}$ [12].

An ideally absorbing detector is used and the source is assumed to have
constant target power (product of tube voltage and current is constant). Its
spectrum is a bremsstrahlung spectrum given by Kramers’ law, which is
sufficient as a rough approximation. Higher tube voltages correspond to a
higher degree of polychromaticity.

To produce plots that are easy to interpret, we again evaluate the
$\mathit{SNR}$ spectrum at a specific spatial frequency $u_{a}$ and assume an
energy-independent MTF. Applications in real world examples must avoid this
kind of simplification.

The monochromatic $\mathit{SNR}$ curve is from the optimal monochromatic x-ray
spectrum (at 6.1 keV for absorption and 7.0 keV for phase contrast) with the
intensity given by the sum over the corresponding bremsstrahlung spectra.555
To get realistic monochromatic intensities, $\frac{1}{3}$ of the sample
thickness is counted as a x-ray filter. The optimal CEW curve uses the signal
strength as the weight. The optimal DEW curve only includes photons below the
optimal threshold, shown as vertical dashed lines in Fig. 4 (middle). The
”$\propto E$ CEW” is the energy integrating indirect detector case ($w\propto
E$) and assumes efficient conversion (see eq. (5), $c^{2}H_{v}^{2}\gg c$).

Both for photoabsorption and phase contrast signals, only the lower energy
x-ray photons have a sufficiently high $\mathit{SDS}$ to contribute positively
to the $\mathit{SNR}$. This effect is weaker for phase contrast due to a
weaker energy dependency. In this example, the $\mathit{SNR}$ curves are
always ordered: $\propto E$ CEW $<$ no weight $\leq$ optimal DEW $<$ optimal
CEW $<$ monochromatic. The differences get larger for broader x-ray spectra.

Direct detectors (”no weight”) here have an intrinsically higher image quality
than indirect detectors (”$\propto E$ CEW”) and this benefit is larger for
higher degrees of polychromaticity. This difference is caused solely by the
different energy weighting, as both detectors are otherwise assumed to be
perfect absorbers and without additional noise.

Figure 5: Simulated (photo-)absorption $\mathit{SNR}(u_{a})$ curves for
different photoabsorber screen thicknesses $t$ (top). Corresponding values for
the weighted absorption average (upper middle) and energy-dependent absorption
efficiency (lower middle). The average detector x-ray absorption for the $t$
with the maximal $\mathit{SNR}$ is shown for different tube voltages (bottom).

It can be seen that weighting down (CEW) or not detecting (DEW) specific
photons can increase $\mathit{SNR}$ by large factors. Note that the
monochromatic $\mathit{SNR}$ is directly proportional to the cumulative
intensity of the source spectrum.

The following simple rules for energy weighting can be seen in Fig. 4: For
optimal CEW, every additional photon increases $\mathit{SNR}$. For optimal DEW
no additional photon decreases $\mathit{SNR}$. For all other weightings,
additional photons can decrease or increase $\mathit{SNR}$.

If we use this example for a $\mathit{SNR}$ optimization, we can see that very
low tube voltages would be optimal if CEW or DEW cannot be implemented.

We could use the optimal CEW case as $\mathit{SNR}_{\text{ideal}}$ in eq. (3)
and the $\propto E$ CEW as $\mathit{SNR}_{\text{detected}}$ to compute an
accurate polychromatic $\mathit{DQE}$ for an indirect detector with ideal
absorption and efficient conversion. It has $\mathit{DQE}\ll 1$ for broad
spectra due to its energy weighting. Computing a polychromatically averaged
$\mathit{DQE}$ using eq. (4) without consideration of the signal strengths
would give $\mathit{DQE}=1$. This is where sample-independent $\mathit{DQE}$
models fail.

In real applications, samples usually have varying thicknesses and a setup
must be optimized to give high $\mathit{SNR}$ at a combination of thicknesses.

###  4.3 Optimal Screen Thickness

$\mathit{SNR}(u_{a})$ curves for different thicknesses $t$ of a x-ray
absorbing screen (e.g. scintillator) are shown in Fig. 5. The screen is
approximated as an idealized photoabsorber with an absorption constant
$\propto E^{-3}$ (no absorption edges). All other properties and the sample
are the same as above and a source spectrum with a tube voltage of 38 kV is
used in the three upper plots. The average absorption values are computed as
the weighted detected spectral average of the x-ray absorption efficiency of
the detection screen. In real cases, the detector MTF depends on the x-ray
detection position within the screen and the intensity distribution within the
screen is energy-dependent.

We can see that for the ”no weight” and the ”$\propto E$ CEW” cases, there is
an optimal screen thickness. Using a thicker screen reduces the $\mathit{SNR}$
due to the lower $\mathit{SDS}$ of the additionally detected higher energy
photons. The ”optimal CEW” and ”optimal DEW” methods however prevent this
effect. In addition, a thinner screen may have a better MTF (this effect is
not simulated here).

In Fig. 5 (bottom) the values of the average x-ray absorption of the detector
for the $t$ with the maximal $\mathit{SNR}$ are shown. It can be seen that in
the absence of optimal CEW or DEW, lower average absorption values are optimal
for higher degrees of polychromaticity. This is an application of the DEW, as
discussed before.

The average absorption for ”no weight” is identical with a polychromatically
averaged $\mathit{DQE}$ (eq. (3)). Thus the ”no weight” and the ”$\propto E$
CEW” cases have an optimal $\mathit{SNR}$ at a specific value of the
polychromatic $\mathit{DQE}$ which is different from the maximal
$\mathit{DQE}$. Raising the $\mathit{DQE}$ usually has other costs (e.g. worse
MTF), so that in this case, increasing the $\mathit{DQE}$ beyond its optimum
can have direct and indirect disadvantages. The ”optimal CEW” and ”optimal
DEW” cases do not benefit significantly if the detector thickness is increased
beyond this point.

##  5 Discussion

### Different Use Cases

The technical capabilities of different x-ray imaging setups vary greatly. The
samples are also very different: Imaging e.g. a
$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick metallic sample at
$100\text{\,}\mathrm{nm}$ resolution or a $500\text{\,}\mathrm{mm}$ organic
sample at $500\text{\,}\mathrm{\SIUnitSymbolMicro m}$ resolution requires very
different setups.

High resolution x-ray imaging (typical sampling $<$
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$) generally uses lower energies
(e.g. $<$ 10 keV) due to the higher x-ray transparency of thinner samples.
Detecting smaller structures also intrinsically requires a higher
$\mathit{SNR}$, so its optimization is generally more important for reducing
measurement time than for coarser structures. The signal in high resolution
x-ray imaging mainly consists of photoabsorption and phase contrast. In low
resolution x-ray imaging, larger samples are investigated. They require higher
x-ray energies. For the latter, Compton scattering produces a significant part
of the absorption signal.

Within the range of x-ray imaging applications, medical imaging represents a
very narrow range. Models for image quality which may be appropriate for these
cases can be inappropriate in a different context, e.g. in material science.
Imaging scientists may therefore need to use models different from those used
in medical imaging (e.g. DQE).

$\mathit{SNR}$ spectra can be used to optimize image quality including phase
contrast effects [12, 14]. The latter can be interpreted as a physical
highpass filter [15]. For x-ray imaging setups with (sub-)micrometer
resolution, phase contrast is often the strongest contribution to the detected
signal.

If one imaging device is used for imaging very different samples, some
compromise must be made on the setup components. Also, if one sample has very
different attenuation lengths along the beam direction, a compromise must be
made to optimize image quality for the corresponding different x-ray spectra.

It is not possible to derive general rules that apply to every case. For any
specific use case, optimizing the $\mathit{SNR}$ spectrum by measurements is
required. For a set of similar use cases, one can derive general rules. Such a
set can e.g. be imaging with a photoabsorption signal and a bremsstrahlung
source. A general rule in this example is the fact that only lower x-ray
energies contribute to image quality while higher energies deteriorate it.
Additionally, simulations are an important method to discover ways in which an
imaging setup may be optimized.

### Structure Size

Using $\mathit{SNR}$ spectra implies that sample structure (of a specific
size) is of interest. Because $\mathit{SNR}(u)$ usually strongly decreases to
higher $u$ (smaller structures), larger structures are always detected well in
a specific imaging measurement. Smaller structures can only be detected down
to a minimal size (”spatial resolution”) which is given by the properties of
the setup and the measurement configuration. For example, the measurement time
can be increased to resolve smaller details.

Depending on the structure size of interest, different imaging setups may be
optimal. This fact is modeled well by $\mathit{SNR}$ spectra. Optimal
detection depends on what is intended to be detected, which in itself is a
decision which needs to be made by the person doing the investigation.

### Temporal Measures

Temporal quantities ($\mathit{SNR}_{t}$ or $\mathit{DQE}_{t}$) are defined as
properties within a single pixel. Due to the fact that the signal in imaging
almost always spans multiple pixels, temporal quantities have no reliable
quantitative relation to image quality. Additionally, the value of a temporal
quantity depends on the arbitrary choice of the area $A$ of one pixel (or
volume and voxel in 3D):

$\mathit{SNR}_{t}\propto\sqrt{A}$ (49)

Detectors with different pixel sizes thus have intrinsically different
$\mathit{SNR}_{t}$, even if the image quality is identical.666This can of
course be easily corrected by normalizing temporal $\mathit{SNR}$ to pixel
area, which is seldom done. This is caused by the fact that the incident
intensity per detector area is independent of the choice of the pixel size and
that temporal $\mathit{SNR}$ is proportional to the square root of the
intensity. $A$ is the area that corresponds to the sampling of the object, not
the sampling of the detector. Changing the x-ray magnification $M$ thus
distorts the temporal $\mathit{SNR}$, which unlike image quality is
independent of x-ray magnification:

$\mathit{SNR}_{t}=\text{const}(M)$ (50)

The actual dependence of the image quality is described by:

$\mathit{SNR}(u)\propto M^{2}$ (51)

which can be confirmed by measuring $\mathit{SNR}(u)$. This is due to the fact
that image quality depends on the detected intensity per object area (not per
detector area), which increases with higher $M$. The smaller area of a object
pixel ($A\propto M^{-2}$) at higher magnifications and the constant intensity
per pixel leads to an increase of the intensity density.

For indirect detectors, temporal measures also neglect the noise correlations
that exist (see section 2.4).

The measurement error in a single pixel does not solely determine image
quality. Instead, image quality depends on the image and what we want to
detect within. The use of temporal image quality measures in an imaging
context must therefore be considered very carefully, because it can easily
result in errors.

### Optimizing Imaging Setups

Designing a good imaging setup or operating it well is conceptually different
in a polychromatic use case than in a monochromatic use case. For
monochromatic imaging, an optimal setup generates the highest possible
intensity and detects all photons. In the polychromatic case one can apply a
CEW and detecting or generating less photons may increase the $\mathit{SNR}$.

If optimal CEW is not possible, then the optimal polychromatic setup needs to
detect as many photons as possible from some photon classes and no photons
from the other photon classes, see eq. (24). Depending on the differentiating
criterion, this may be very difficult to achieve. While a detector can be
designed to less efficiently detect x-ray photons of a specific energy range,
this is much more difficult for x-ray photons that were absorbed at a position
in the screen where the MTF is worse (e.g. higher distance from the focus
plane).

Counting detectors which can set an upper limit to the x-ray energy of the
counted detection events are ideally suited for DEW. Setting the upper limit
in such a way that photons with comparatively low $\mathit{SDS}$ are not
counted may potentially result in a large increase in $\mathit{SNR}$ (see Fig.
4).

For an indirect detector, the $\mathit{SNR}$ can be increased by choosing a
thinner screen for lower energy imaging – even if the absorption efficiency at
the relevant energy range may then be slightly lower (see Fig. 5). This effect
is separate from the better MTF of the thinner screen. In many imaging setups,
the screen is already thin enough for its use cases due to technical
limitations or MTF considerations.

### Testable Predictions

An image quality optimization based on $\mathit{SNR}$ spectra or temporal
$\mathit{SNR}$ makes predictions which can be tested directly. This is done
with a direct measurement, see A or [14] for a method. A $\mathit{SNR}$
optimization should ideally be based on measurements for the specific use
case. Using simulations instead requires that the simulation method was
thoroughly tested with direct measurements on the same or a sufficiently
similar device. This makes it possible to notice errors in the theory or in
the assumptions made, which is a basic requirement if one wants to rely on
such a model [16].

Concerning $\mathit{DQE}$, note that optimizing $\mathit{SNR}$ may lead to a
different optimal setting than optimizing polychromatically averaged
$\mathit{DQE}$. This difference stems from the fact that in practice
$\mathit{DQE}$ measurements must use assumptions that may be wrong.

A model can both be tested by physical measurements and also by checking if it
is without internal contradictions. Here, internal contradictions mainly mean
that the model does not actually describe what it is supposed to describe.
This is for example the case when temporal $\mathit{SNR}$ is used in imaging.

##  6 Conclusions

The concept of signal detection allows a systematic investigation of physical
effects that result from varying signal or detection properties for the
different detected photons. Image quality measures like the $\mathit{DQE}(u)$
falsely assume these differences to be unimportant and are therefore unable to
model cases where they are.

Temporal measures like the $\mathit{SNR}_{t}$ should principally not be used
in imaging due to the fact that an image consists of more than one pixel.

We use $\mathit{SNR}$ spectra as a quantitative model for image quality within
the framework of signal detection. This allows a robust optimization of image
quality in polychromatic x-ray imaging. The predictions of this model can and
should be tested with direct measurements. An image quality optimization on an
existing imaging device can be done with direct measurements without the need
to correctly model its physical properties.

The main disadvantage $\mathit{SNR}$ spectra have is that they are not an
absolute quantitative measure. The actual value depends on the object. While
this is necessary for an accurate model, a quantity that depends only on the
imaging setup is more convenient. This is solved by simplifying the
$\mathit{SNR}$ spectra to a quantity like the detection effectiveness
$\mathit{DE}(u)$, eq. (18), at the cost of a more complex measurement process.
The $\mathit{DE}(u)$ is not as convenient as an absolutely normalized quantity
like the $\mathit{DQE}(u)$, but more robust and accurate.

For polychromatic imaging, we have shown that in a case with varying signal
strengths, reducing the detected intensity or weight of photons with a
relatively low signal detection strength ($\approx$ $\mathit{SDS}$, eq. (26))
can increase the $\mathit{SNR}$ and therefore the image quality. While the
absolute signal strength does not influence an image quality optimization,
relative strength does. Of course, increasing the detected intensity of
photons with a high $\mathit{SDS}$ does increases image quality. In principle,
this result can be applied to any imaging setup.

The model developed here is especially useful for high resolution x-ray
imaging and phase contrast imaging. We found that depending on the
application, direct detectors can have an intrinsically higher $\mathit{SNR}$
than indirect detectors, due to the different energy weightings.

It is likely that much work analyzing image quality using the $\mathit{DQE}$
or temporal $\mathit{SNR}$ will transfer to the concept of signal detection
without large differences in the conclusion. In these cases, the benefit will
mainly be that the underlying assumptions become clear, which in turn may
prevent false generalizations. In other cases, using the signal detection
model is required to accurately model and optimize the image quality of a
polychromatic x-ray imaging device.

### Acknowledgments

We would like to thank Theobald Fuchs for test reading and feedback and Alison
Haydock for language help.

### Supplementary Material

  * $\circ$

The three Jupyter notebooks (Python) that were used to generate figures,
including the simulations.

  * $\circ$

Python code for a $\mathit{SNR}$ spectra evaluation from measurements as
described in the A. Includes examples in two Jupyter notebooks.

## References

  * Yaffe and Rowlands [1997] MJ Yaffe and JA Rowlands.  X-ray detectors for digital radiography.  _Phys Med Biol_ , 42(1):1, 1997.  10.1088/0031-9155/42/1/001.
  * Cunningham and Shaw [1999] Ian A Cunningham and Rodney Shaw.  Signal-to-noise optimization of medical imaging systems.  _J Opt Soc Am A_ , 16(3):621–632, 1999.  10.1364/josaa.16.000621.
  * Stampanoni et al. [2002] Marco Stampanoni, Gunther Borchert, Peter Wyss, Rafael Abela, Bruce Patterson, Steven Hunt, Detlef Vermeulen, and Peter Rüegsegger.  High resolution x-ray detector for synchrotron-based microtomography.  _Nucl Intrum Meth A_ , 491(1):291–301, 2002\.  10.1016/S0168-9002(02)01167.
  * Siewerdsen et al. [2002] J. H. Siewerdsen, I. A. Cunningham, and D. A. Jaffray.  A framework for noise-power spectrum analysis of multidimensional images.  _Med Phys_ , 29(11):2655–2671, 10 2002.  10.1118/1.1513158.
  * Tanguay et al. [2015] Jesse Tanguay, Seungman Yun, Ho Kyung Kim, and Ian A Cunningham.  Detective quantum efficiency of photon-counting x-ray detectors.  _Med Phys_ , 42(1):491–509, 2015.  10.1118/1.4903503.
  * Siewerdsen et al. [2004] J. H. Siewerdsen, A. M. Waese, D. J. Moseley, S. Richard, and D. A. Jaffray.  Spektr: A computational tool for x-ray spectral analysis and imaging system optimization.  _Med Phys_ , 31(11):3057–3067, 10 2004.  10.1118/1.1758350.
  * Illers et al. [2005] Hartmut Illers, Egbert Buhr, and Christoph Hoeschen.  Measurement of the detective quantum efficiency (dqe) of digital x-ray detectors according to the novel standard iec 62220-1.  _Radiat Prot Dosimetry_ , 114(1-3):39–44, 2005\.  10.1093/rpd/nch507.
  * IEC [2003] Iec 62220-1: Medical electrical equipment – characteristics of digital x-ray imaging devices – part 1: Determination of the detective quantum efficiency.  IEC, 10 2003.
  * Jähne [2010] Bernd Jähne.  Emva 1288 standard for machine vision: Objective specification of vital camera data.  _Optik & Photonik_, 5(1):53–54, 2010.  10.1002/opph.201190082.
  * Cunningham et al. [1994] Ian Alexander Cunningham, MS Westmore, and A Fenster.  A spatial-frequency dependent quantum accounting diagram and detective quantum efficiency model of signal and noise propagation in cascaded imaging systems.  _Med Phys_ , 21(3):417–427, 1994.  10.1118/1.597401.
  * Bracewell [1986] R.N. Bracewell.  _The Fourier Transform and Its Applications_.  McGraw-Hill, 1986.  ISBN 0-306-48187-1.
  * Paganin [2006] D. M. Paganin.  _Coherent X-Ray Optics_.  Oxford University Press, 2006.  ISBN 9780199673865.
  * Tapiovaara and Wagner [1985] M J Tapiovaara and R Wagner.  SNR and DQE analysis of broad spectrum x-ray imaging.  _Phys Med Biol_ , 30(6):519–529, 6 1985.  10.1088/0031-9155/30/6/002.
  * Ullherr et al. [2018] M Ullherr, A Balles, C Fella, and S Zabler.  Using measurements of the spatial snr to optimize phase contrast x-ray imaging.  _Nucl Intrum Meth A_ , 877:44–50, 2018.  10.1016/j.nima.2017.09.044.
  * Paganin et al. [2002] D. Paganin, SC Mayo, TE Gureyev, PR Miller, and SW Wilkins.  Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object.  _J Microsc_ , 206(1):33–40, 2002.  10.1046/j.1365-2818.2002.01010.x.
  * Popper [2005] Karl Popper.  _The logic of scientific discovery_.  Routledge, 2005.  10.4324/9780203994627.

##  Appendix A Measuring SNR Spectra

An important aspect of the $\mathit{SNR}$ Spectra is that they can directly be
measured. The method described in [14] can be used for this purpose—also in
applications other than phase contrast. In the following, we will give a short
description of how to apply this method.

The measurement consists of using an appropriate test object of specific
material composition and thickness as sample and acquiring a series of $K$
images $\\{d_{k}\\}_{k}$ with this sample. The images can either be
projections or CT reconstructions. Individual images must differ only in their
noise realizations. The individual exposure times can be shorter than those
used in practice; the images should be visibly noisy. The averaged image is
defined as:

$d_{\text{avg}}(x)=\frac{1}{K}\sum_{k=1}^{K}d_{k}(x)$ (52)

If $D(u)$ is the average power spectrum of the $d_{k}$, $D_{\text{avg}}(u)$
the power spectrum of $d_{\text{avg}}(x)$, $\tau$ the integration time of the
x-ray detector, then the $\mathit{SNR}$ spectrum for this setup can be
calculated as:

$\mathit{SNR}_{\tau}(u)=\frac{D_{\text{avg}}(u)-\tfrac{1}{K}D(u)-(1-\tfrac{1}{K})A(u)}{D(u)-D_{\text{avg}}(u)}$
(53)

where $A(u)$ is the power spectrum of the known image artifacts (e.g.
reference image noise). This formula works for $K\geq 2$ but gives better
results for $K>$ 20\. Note that for computing the power spectra from real
images with a FFT, using an appropriate window function is necessary to get an
accurate result.

If a reference image is used to normalize the intensity of the projection
image, the noise power spectrum of a single reference image can be measured
similarly as above with:

$\mathit{N}_{\tau,\text{ref}}(u)=\frac{D_{\text{ref}}(u)-D_{\text{avg,\,ref}}(u)}{1-\tfrac{1}{K_{\text{ref}}}}$
(54)

This is similar to the denominator of eq. (53) except for the factor which
canceled out. If the average of the ref images is used for the normalization,
it gives a contribution to $A(u)$ of the form:

$A(u)=...+\frac{\mathit{N}_{\tau,\text{ref}}(u)}{K_{\text{ref}}}$ (55)

Detector imperfections or other image artifacts not modeled by $A(u)$ can
falsely increase the calculated signal and $\mathit{SNR}$. These signal
artifacts usually appear as a lower limit to $\mathit{SNR}$ at the higher
spatial frequencies where the actual $\mathit{SNR}$ becomes very small. For a
good measurement, this lower limit can be as low as 10-3, which is not a
practical problem.

$\mathit{SNR}$ spectra measurements from different samples (”test phantoms”)
can generally not be compared and their absolute values have limited
usefulness. The reason is that different amounts of sample structure lead to
different $\mathit{SNR}$ spectra. In many cases, comparability is not needed
or the results can be stated as a relative difference between different
conditions for the same sample.

Otherwise, the detection effectiveness spectrum defined in eq. (18) can be
used. This quantity can be determined for an experimental setup if the
$\mathit{SNR}$ spectrum is measured and the object spectrum $P(u)$ is computed
for the object used. For full comparability, the sample thicknesses must be
approximately the same for all measurements that need to be compared. It is
possible to choose a test phantom for which $P(u)$ can easily be computed,
e.g. one consisting of balls. This implies that the $\mathit{SNR}$ measurement
is done for the CT reconstruction image.

Test phantoms may need to be combined with x-ray filters of different
thickness to simulate different sample thicknesses. The cumulative thickness
of all objects in the beam is an integral part of the definition of the test
phantom. Any test phantom should satisfy the following properties:

  * $\circ$

Its power spectrum should not change with small differences ($<$
$1\text{\,}\mathrm{\SIUnitSymbolDegree}$ or $<$ cone angle) in the orientation
of the sample to the beam. Round shapes are ideal, e.g. balls with a diameter
of 20-50 voxels.

  * $\circ$

It should represent a type of sample or a category of samples which will
realistically be used.

Another possibility of obtaining comparable measurements is to define
standardized test phantoms. Examples of such definitions are: (1) A
geometrically defined object of a specific material. (2) A random object with
a defined statistic, e.g. a size distribution of balls and a specific
container geometry.

Note that robust measurements of the temporal $\mathit{SNR}$ require the same
type of measured data but have a slightly different evaluation. Estimating
temporal $\mathit{SNR}$ from a region with constant gray values is prone to
errors.


# Hamiltonian mechanics of “magnetic” solitons in two-component Bose-Einstein
condensates

A. M. Kamchatnov  Institute of Spectroscopy, Russian Academy of Sciences,
Troitsk, Moscow, 108840, Russia  Skolkovo Institute of Science and Technology,
Skolkovo, Moscow, 143026, Russia

###### Abstract

We consider motion of a “magnetic” soliton in two-component condensates along
a non-uniform and time-dependent backgrounds in framework of the Hamiltonian
mechanics. Our approach is based on generalization of Stokes’ remark that
soliton’s velocity is related with its inverse half-width by the dispersion
law for linear waves continued to the region of complex wave numbers. We
obtain expressions for the canonical momentum and the Hamiltonian as functions
of soliton’s velocity and transform the Hamilton equations to the Newton-like
equation. The theory is illustrated by several examples of concrete soliton’s
dynamics.

Dedicated to the memory of David Kaup

##  I Introduction

The idea that solitons behave under action of smooth enough external forces
similar to point-like particles of classical mechanics is well known. In
particular, a number of perturbation methods is based on this qualitative
picture [1, 2, 3, 4] (see also Refs. [5, 6] and references therein) and these
theories are confirmed by numerical and real experiments. The standard
perturbation theory is quite involved in case of dark solitons propagating
along a non-uniform background [7, 8, 9] when separation of the soliton’s
dynamics from dynamics of the surrounding background is not obvious: a moving
soliton generates a ‘shelf’ in the background and this leads to an essential
back reaction on the soliton’s motion. For example, the frequency of
oscillations of a dark soliton in a Bose-Einstein condensate (BEC) confined in
a harmonic trap with the frequency $\omega_{0}$ is equal to
$\omega_{0}/\sqrt{2}$ whereas the center of gravity of BEC still oscillates
with the trap frequency [8]. This means that the soliton’s motion generates a
counterflow in the surrounding condensate which results in an essential
redistribution of its mass. In this sense, motion of dark solitons can be
likened [10] to motion of quasi-particles in solid-state physics where
dynamical properties of quasi-particles are determined by their interaction
with the surrounding matter. As a result, derivation of equations of motion of
such solitonic quasi-particles by the standard perturbation theory is,
generally speaking, quite a difficult problem.

A simpler method has been recently suggested in Ref. [11]. It is based on the
assumption that if the nonlinear wave equation is Hamiltonian, then its
reduction to the case of the soliton’s motion along a smooth time-dependent
background remains also Hamiltonian. The condition that the soliton’s dynamics
is Hamiltonian, combined with old Stokes’ remark [12] that the soliton’s
velocity is related with its inverse half-width by the dispersion relation for
linear waves, allows one to find expressions for the canonical momentum and
the Hamiltonian of the point-like soliton. The method was applied to the KdV
[11], NLS [13, 14], DNLS [15] equations and analytical results were confirmed
by numerical solutions for concrete situations.

In this paper, we will obtain the equations of motion for so-called “magnetic”
solitons in a two-component BEC. Such soliton solutions were obtained in Ref.
[16] in the limit of small difference between the inter- and intra-components
interaction. In this case, the BEC components move in opposite directions
without transfer of the total density, so the dynamics is completely
determined by the counterflow of the BEC components. We obtain the expressions
for the canonical momentum and the Hamiltonian of the soliton dynamics and
transform the Hamilton equations to the Newton-like equation for soliton’s
motion. Since the dynamics is determined by the gradients of the background
flow, we assume that the Newton equation remains valid when these gradients
are resulted from the background dynamics with account of external potentials.
The theory is illustrated by concrete examples of the magnetic soliton
dynamics.

##  II Basic equations

A two-component condensate is described by the system of Gross-Pitaevskii
equations

$\begin{split}&i\psi_{1,t}+\frac{1}{2}\psi_{1,xx}-(\alpha_{1}|\psi_{1}|^{2}+\alpha_{2}|\psi_{2}|^{2})\psi_{1}=V_{1}\psi_{1},\\\
&i\psi_{2,t}+\frac{1}{2}\psi_{2,xx}-(\alpha_{2}|\psi_{1}|^{2}+\alpha_{1}|\psi_{2}|^{2})\psi_{2}=V_{2}\psi_{2},\end{split}$
(1)

where $\psi_{1,2}$ are wave functions of two components, $\alpha_{1}$ is the
constant of interaction between atoms in each component, and $\alpha_{2}$ is
the constant of interaction between atoms which belong to different
components. We assume that

$\alpha_{1}>\alpha_{2}>0,$ (2)

what corresponds to the repulsive interaction between atoms and the BEC is a
mixture of the components which do not separate from each other in space (see,
e.g., [17]). $V_{1}$ and $V_{2}$ are the external (trap) potentials acting on
the first and the second component, correspondingly.

It is convenient to introduce pseudospin variables (see, e.g., [18])

$\left(\begin{array}[]{c}\psi_{1}\\\ \psi_{2}\\\
\end{array}\right)=\sqrt{\rho}\,e^{i\Phi/2}\left(\begin{array}[]{c}\cos\frac{\theta}{2}\,e^{-i\phi/2}\,e^{-i\mu_{1}t}\\\
\sin\frac{\theta}{2}\,e^{i\phi/2}\,e^{-i\mu_{2}t}\\\ \end{array}\right),$ (3)

where $\rho_{1}=\rho\cos^{2}(\theta/2),\rho_{2}=\rho\sin^{2}(\theta/2)$ are
the densities of the components, so that $\rho=\rho_{1}+\rho_{2}$ is the total
density of BEC. The functions $\Phi$ and $\phi$ play the roles of potentials
of velocities of the mutual and relative motions of the components;
consequently the components’ velocities are equal to

$v_{1}=\frac{1}{2}(\Phi_{x}-\phi_{x}),\qquad
v_{2}=\frac{1}{2}(\Phi_{x}+\phi_{x}).$ (4)

In the Thomas-Fermi approximations, Eqs. (1) give the stationary distributions
of the densities of the components confined in a trap with the potentials
$V_{1}$ and $V_{2}$,

$\begin{split}\rho_{1}=\rho_{1,0}(x)=\frac{\alpha_{1}(\mu_{1}-V_{1}(x))-\alpha_{2}(\mu_{2}-V_{2}(x))}{\alpha_{1}^{2}-\alpha_{2}^{2}}),\\\
\rho_{2}=\rho_{2,0}(x)=\frac{\alpha_{1}(\mu_{2}-V_{2}(x))-\alpha_{2}(\mu_{1}-V_{1}(x))}{\alpha_{1}^{2}-\alpha_{2}^{2}};\end{split}$
(5)

then the total density is equal to

$\rho=\rho_{0}(x)=\frac{1}{\alpha_{1}+\alpha_{2}}(\mu-V_{1}(x)-V_{2}(x)),$ (6)

where $\mu=\mu_{1}+\mu_{2}$.

Substitution of Eqs. (3) into Gross-Pitaevskii Equations (1) yields [19]

$\begin{split}&\rho_{t}+\frac{1}{2}\left[\rho(\Phi_{x}-\cos\theta\phi_{x})\right]_{x}=0,\\\
&\Phi_{t}+\frac{\rho_{x}^{2}}{4\rho^{2}}-\frac{\rho_{xx}}{2\rho}-\frac{\cot\theta}{2\rho}(\rho\,\theta_{x})_{x}+\frac{1}{4}(\Phi_{x}^{2}+\phi_{x}^{2}+\theta_{x}^{2})+(\alpha_{1}+\alpha_{2})\rho=-(V_{1}+V_{2}),\\\
&\theta_{t}+\frac{1}{2\rho}(\phi_{x}\rho\sin\theta)_{x}+\frac{1}{2}\Phi_{x}\theta_{x}=0,\\\
&\phi_{t}-\frac{1}{2\rho\sin\theta}(\rho\,\theta_{x})_{x}+\frac{1}{2}\Phi_{x}\phi_{x}-(\alpha_{1}-\alpha_{2})\rho\cos\theta=V_{1}-V_{2}.\end{split}$
(7)

We are interested in dynamics of “magnetic” solitons of relative motion of the
components [20] and this solitons correspond to fast changes of the relative
density (that is of the $\theta$-variable) and of the relative velocity
$v=v_{2}-v_{1}=\phi_{x}$. The characteristic width $\kappa$ of such a soliton
can be easily estimated from the last equation,
($(\rho_{0}\kappa^{2})^{-1}\sim(\alpha_{1}-\alpha_{2})\rho$), that is in our
non-dimensional variables

$\kappa\sim(\rho_{0}\delta g)^{-1/2},$ (8)

where we introduced the notation $\delta g=\alpha_{1}-\alpha_{2}$. Of course,
it is supposed that the soliton’s width is much smaller than the
characteristic size of distributions (5) and (6). The second equation (7) can
be written with account of Eq. (6) in the form

$\Big{\\{}\frac{\rho_{X}^{2}}{4\rho^{2}}-\frac{\rho_{XX}}{2\rho}-\frac{\cot\theta}{2\rho}(\rho\,\theta_{X})_{X}+\frac{1}{4}(\Phi_{X}^{2}+\phi_{X}^{2}+\theta_{X}^{2})\Big{\\}}\kappa^{-2}=2g(\rho_{0}-\rho),$

where $2g=\alpha_{1}+\alpha_{2}$ and $X=x/\kappa$, that is $X\sim 1$ within
the soliton’s width. The expression inside the curly brackets is of order of
magnitude of unity and it is multiplied by the factor $\kappa^{-2}\sim\delta
g$ which is very small if the difference $\delta g$ of the interaction
constants is small. Consequently, in this limit discussed in Ref. [16], the
total density can be considered as constant within the soliton’s width: there
is no transfer of total mass due to propagation of a magnetic soliton, i.e.
$\rho=\rho_{0}$. Thus, we arrive at the picture of point-like soliton
propagating along smooth distributions (5) and (6).

If the potentials are stationary, then $\rho_{t}=0$, so the first Eq. (7)
gives $\rho(\Phi_{x}-\cos\theta\phi_{x})=j=\mathrm{const}$. We will confine
ourselves to situations when there is no flow of the total density, $j=0$, so

$\Phi_{x}=\cos\theta\phi_{x},$ (9)

and the remaining Eqs. (7) reduce to

$\begin{split}&\theta_{t}+\theta_{x}\phi_{x}\cos\theta+\frac{1}{2}\sin\theta\phi_{xx}=-\frac{\rho_{0,x}}{2\rho_{0}}{\phi_{x}}{\sin\theta},\\\
&\phi_{t}-\delta
g\,\rho_{0}\cos\theta+\frac{1}{2}\cos\theta\phi_{x}^{2}-\frac{\theta_{xx}}{2\sin\theta}=V_{1}-V_{2}+\frac{\rho_{0,x}}{2\rho_{0}}\frac{\theta_{x}}{\sin\theta}.\end{split}$
(10)

Here $\theta$ and $\phi$ change at distances of order of magnitude of
soliton’s width, i.e., $\rho_{0,x}/\rho_{0}\ll\theta_{x},\phi_{x}$ and terms
with $\rho_{0,x}/\rho_{0}$ can be neglected. Besides that, to simplify
notation and comparison with other publications, it is conveniens to introduce
a new space coordinate $x^{\prime}=\sqrt{2}x$. Then we get

$\begin{split}&\theta_{t}+2\theta_{x}\phi_{x}\cos\theta+\sin\theta\phi_{xx}=0,\\\
&\phi_{t}-(\delta
g\,\rho_{0}-\phi_{x}^{2})\cos\theta-\frac{\theta_{xx}}{\sin\theta}=V_{1}-V_{2},\end{split}$
(11)

where we still write $x$ instead of $x^{\prime}$. The potentials
$V_{1},V_{2},$ and the total density $\rho_{0}$ are slow functions of $x$.
These equations can be obtained from the Hamilton principle with the
Lagrangian density

$\mathcal{L}=\cos\theta\phi_{t}-\frac{1}{2}\left[\theta_{x}^{2}-\left(\delta
g\,\rho_{0}-\phi_{x}^{2}\right)\sin^{2}\theta\right]-(V_{1}-V_{2})\cos\theta,$
(12)

where $V_{1},V_{2}$, and $\rho_{0}$ are assumed constant. It coincides with
the Lagrangian of a magnetic material with an easy-plane anisotropy (see,
e.g., [21]). If we denote

$\cos\theta=w,\quad\phi_{x}=v,\quad\widetilde{\rho}=2\delta g\,\rho_{0}$ (13)

then we cast Eqs. (11) to the form [19, 22]

$\begin{split}&w_{t}-\left[v(1-w^{2})\right]_{x}=0,\\\
&v_{t}-\left[w(\widetilde{\rho}-v^{2})\right]_{x}+\left[\frac{1}{\sqrt{1-w^{2}}}\left(\frac{w_{x}}{\sqrt{1-w^{2}}}\right)_{x}\right]_{x}=V_{1,x}-V_{2,x}.\end{split}$
(14)

These are nonlinear equations for polarization waves propagating along the
non-uniform stationary background with the total density distribution (6);
Eqs. (5) determine the smooth background distribution of the variable $w$. Our
task is to reduce Eqs. (14) to equations of motion of a point-like magnetic
soliton.

##  III Dispersion relation of linear waves, dispersionless limit and soliton
solution

Let a small-amplitude harmonic wave propagate along a uniform condensate with
constant values of $w=w_{0},v=v_{0}$ and zero external potentials
$V_{1}=V_{2}=0$. We linearize Eqs. (14) with respect to small deviations
$w^{\prime}=w-w_{0},v^{\prime}=v-v_{0}$ from this uniform state and obtain

$\begin{split}&w^{\prime}_{t}+2v_{0}w_{0}w^{\prime}_{x}-(1-w_{0}^{2})v^{\prime}_{x}=0,\\\
&v^{\prime}_{t}-(\widetilde{\rho}-v_{0}^{2})w^{\prime}_{x}+\frac{w^{\prime}_{xxx}}{1-w_{0}^{2}}+2v_{0}w_{0}v^{\prime}_{x}=0.\end{split}$
(15)

This system readily gives the dispersion relation for harmonic waves
$w^{\prime},v^{\prime}\propto\exp[i(kx-\omega t)]$:

$\omega(k)=k\left(2v_{0}w_{0}\pm\sqrt{(\widetilde{\rho}-v_{0}^{2})(1-w_{0}^{2})+k^{2}}\right).$
(16)

By definition we have $|w_{0}|<1$, so the uniform state is modulationally
unstable for large enough values of the relative velocity,
$|v_{0}|>\sqrt{\widetilde{\rho}}$, of the components. We will confine
ourselves to the modulationally stable situations with
$|v_{0}|<\sqrt{\widetilde{\rho}}$.

Even for zero potentials $V_{1}=V_{2}=0$ the background state can be non-
uniform and time-dependent, if a large-scale wave propagates through BEC, so
that the wave variables are give by the functions $w=w(x,t)$ and $v=v(x,t)$.
If the characteristic size of these distributions is much greater than the
soliton’s width $\sim\kappa^{-1}$, then we can neglect the dispersion term in
the second equation (14) and arrive at the equations of the dispersionless
(hydrodynamic) approximation

$w_{t}-[v(1-w^{2})]_{x}=0,\qquad v_{t}-[w(\widetilde{\rho}-v^{2})]_{x}=0.$
(17)

These are Ovsyannikov equations [23] derived long ago for internal waves in
two-fluid systems. Characteristic velocities of this system are equal to

$v_{\pm}=2vw\pm\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})}.$ (18)

Eqs. (17) can be cast to the diagonal Riemann form

$\frac{\partial r_{+}}{\partial t}+v_{+}\frac{\partial r_{+}}{\partial
x}=0,\qquad\frac{\partial r_{-}}{\partial t}+v_{-}\frac{\partial
r_{-}}{\partial x}=0$ (19)

for the Riemann invariants

$r_{\pm}=vw\pm\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})}.$ (20)

The characteristic velocities are expressed in terms of the Riemann invariants
by simple formulas

$v_{+}=\frac{1}{2}(3r_{+}+r_{-}),\qquad v_{-}=\frac{1}{2}(r_{+}+3r_{-}).$ (21)

The simplest soliton solution of Eqs. (14) (with $V_{1}=V_{2}=0$) was found in
Ref. [16] and it corresponds to zero background relative velocity, $v_{0}=0$,
and equal background densities of the components, $w_{0}=0$, far from the
soliton. The general soliton solution with arbitrary values of the background
parameters $v_{0},w_{0}$ was obtained in Ref. [24]. We don’t need here the
explicit form of this solution. Some necessary most important for us relations
can be found without knowledge of the full solution. For example, according to
the old remark of Stokes [12], the exponentially small soliton’s tails
$\propto\exp[\pm\kappa(x-Vt)]$, $V$ being the soliton’s velocity, and the
small-amplitude harmonic waves $\propto\exp[i(kx-\omega t)]$ obey the same
linearized system (15). Consequently, the replacement $k\to i\kappa$
transforms the harmonic wave to the soliton’s tail and the phase velocity
$\omega/k$ of a harmonic wave transforms to the soliton’s velocity,

$V=\frac{\omega(i\kappa)}{i\kappa}.$ (22)

As a result, we obtain the expression for the soliton’s velocity

$V=2vw\pm\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})-\kappa^{2}},$ (23)

where $v,w$ are the local values of the large-scale wave variables obeying the
dispersionless equations (17) and $\kappa$ is the soliton’s inverse half-
width.

Another relationship

$\kappa^{2}=(\widetilde{\rho}-v^{2})(1-w^{2})-(q-vw)^{2}$ (24)

($q$ is a constant parameter) can be obtained using the fact that in
completely integrable system (14) (with $V_{1}=V_{2}=0$) the carrier wave
number of a wave packet propagating along a large-scale wave is only a
function of the local values of $v$ and $w$ given by the formula (see [25])

$k^{2}=(q-vw)^{2}-(\widetilde{\rho}-v^{2})(1-w^{2}),$ (25)

where $q$ is an integration constant determined by the values of $k$ at some
initial point. Again the Stokes transformation yields the formula (24).
Combining Eqs. (24) and (25), we obtain (for the choice of the upper sign and
supposing $q-vw>0$)

$V=\frac{dx}{dt}=vw+q,$ (26)

where $x=x(t)$ is the soliton’s path.

It is instructive to derive this relation by another method based on the
analysis of the Whitham modulation equations at their soliton limit. It was
used in Refs. [26, 27] for consideration of soliton-mean flow interaction and
we generalize here this approach. Periodic solutions of the system (14) with
$V_{1}=V_{2}=0$ are parameterized by four constant parameters $r_{i},$
$i=1,2,3,4,$ which in a modulated wave become slow functions of $x$ and $t$
(see [24]). Their slow evolution is governed by the Whitham modulation
equations

$\frac{\partial r_{i}}{\partial t}+v_{i}\frac{\partial r_{i}}{\partial
x}=0,\quad i=1,2,3,4.$ (27)

At the soliton limit we have $r_{2}=r_{3}$ and these equations reduce to

$\frac{\partial r_{1}}{\partial t}+\frac{1}{2}(3r_{1}+r_{4})\frac{\partial
r_{1}}{\partial x}=0,\qquad\frac{\partial r_{2}}{\partial
t}+\frac{1}{2}(r_{1}+2r_{2}+r_{4})\frac{\partial r_{2}}{\partial
x}=0,\qquad\frac{\partial r_{4}}{\partial
t}+\frac{1}{2}(r_{1}+3r_{4})\frac{\partial r_{4}}{\partial x}=0,$ (28)

and the phase velocity of the periodic solution becomes the soliton’s velocity
equal to

$V=\frac{1}{2}(r_{1}+2r_{2}+r_{4}).$ (29)

Consequently, the soliton’s velocity changes along its path according to the
equation

$\frac{dV}{dt}=\frac{\partial V}{\partial t}+V\frac{\partial V}{\partial
x}=\frac{1}{2}\left[(r_{2}-r_{1})\frac{\partial r_{1}}{\partial
x}-(r_{4}-r_{2})\frac{\partial r_{4}}{\partial x}\right].$ (30)

The first and third Eqs. (28) coincide with Eqs. (19), that is they describe
evolution of the large-scale wave at the soliton’s location and we can
identify $r_{1}=r_{-},r_{4}=r_{+}$, where $r_{\pm}$ are defined bt Eqs. (20).
Substitution of these formulas and of $r_{2}=V-(r_{1}+r_{4})/2=V-vw$ into Eq.
(30) yields

$\frac{dV}{dt}=(1-3w^{2})vv_{x}+(\widetilde{\rho}-3v^{2})vv_{x}+V(vw)_{x}=(vw)_{t}+V(vw)_{x}=\frac{d(vw)}{dt}.$
(31)

Integration of this equation gives at once Eq. (26).

##  IV Hamilton equations

In derivation of the canonical momentum and Hamiltonian, we still assume that
the background wave obeys the dispersionless equations (17) with
$\widetilde{\rho}=\mathrm{const}$. It is convenient to introduce a new
variable $\varphi$ instead of $\kappa$,

$\kappa=\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})}\sin\varphi.$ (32)

Then Eq. (23) takes the form

$V=2vw+\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})}\cos\varphi$ (33)

and substitution of this expression into Eq. (26) gives

$vw+\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})}\cos\varphi=q.$ (34)

We will need a useful formula which follows from this identity. To obtain it,
we differentiate (34) along the soliton’s path $x=x(t)$, so that
$d/dt=\partial/\partial t+V\partial/\partial x$, and eliminate $\partial
v/\partial t$, $\partial w/\partial t$ with the use of Eqs. (17). Then evident
simplifications yield the necessary formula

$\frac{d\varphi}{dt}=-\left[\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})}\right]_{x}\sin\varphi.$
(35)

Now we interpret Eq. (33) as the Hamilton equation

$\frac{dx}{dt}=2vw+\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})}\cos\varphi=\left(\frac{\partial
H}{\partial p}\right)_{x}$ (36)

and assume that the momentum has the form

$p=(\widetilde{\rho}-v^{2})(1-w^{2})f(\varphi),$ (37)

where $f(\varphi)$ is an unknown function. Integration of Eq. (36) gives

$H=2vwp+\left[(\widetilde{\rho}-v^{2})(1-w^{2})\right]^{3/2}\int\cos\varphi
f^{\prime}(\varphi)d\varphi.$ (38)

To find $f(\varphi)$, we turn to another Hamilton equation

$\frac{dp}{dt}=-\left(\frac{\partial H}{\partial x}\right)_{p}.$ (39)

Substitution of Eqs. (37) and (38) gives after evident simplifications

$\begin{split}&\left[(\widetilde{\rho}-v^{2})(1-w^{2})\right]_{x}\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})}f\cos\varphi+(\widetilde{\rho}-v^{2})(1-w^{2})f^{\prime}\frac{d\varphi}{dt}\\\
&=-\frac{3}{2}\sqrt{(\widetilde{\rho}-v^{2})(1-w^{2})}\left[(\widetilde{\rho}-v^{2})(1-w^{2})\right]_{x}\int\cos\varphi
f^{\prime}d\varphi-\left[(\widetilde{\rho}-v^{2})(1-w^{2})\right]^{3/2}\cos\varphi
f^{\prime}\left(\frac{\partial\varphi}{\partial x}\right)_{p}.\end{split}$
(40)

Differentiation of Eq. (37) with respect to $x$ for constant $p$ gives

$\left(\frac{\partial\varphi}{\partial
x}\right)_{p}=-\frac{f\left[(\widetilde{\rho}-v^{2})(1-w^{2})\right]_{x}}{f^{\prime}(\widetilde{\rho}-v^{2})(1-w^{2})}.$
(41)

We substitute Eqs. (35) and (41) into Eq. (40) and arrive at the equation

$f^{\prime}\sin\varphi=3\int^{\varphi}\cos\varphi\,f^{\prime}d\varphi$ (42)

which can be solved at once to give

$f=C(\varphi-\sin\varphi\cos\varphi),$ (43)

where $C$ is an integration constant. Its arbitrariness corresponds to the
invariance of the Hamilton equations with respect to multiplication of $H$ and
$p$ by the same constant factor. We choose $C=2$ to simplify some further
formulas. Then Eqs. (37) and (38) yield after exclusion of $\varphi$ with help
of Eq. (33) the following expressions for the canonical momentum

$p=2(1-w^{2})(\widetilde{\rho}-v^{2})\arccos\frac{\dot{x}-2vw}{\sqrt{(1-w^{2})(\widetilde{\rho}-v^{2})}}-2(\dot{x}-2vw)\sqrt{(1-w^{2})(\widetilde{\rho}-v^{2})-(\dot{x}-2vw)^{2}},$
(44)

and for the Hamiltonian

$H=2vwp+\frac{4}{3}\left[(1-w^{2})(\widetilde{\rho}-v^{2})-(\dot{x}-2vw)^{2}\right]^{3/2},$
(45)

where $x=x(t)$ is the soliton’s trajectory.

It is worth noticing that these canonical momentum and Hamiltonian are related
with those functions in Hamiltonian mechanics of NLS equation solitons [13,
14] by simple formulas. Indeed, the ‘shallow water equations’

$n_{t}+(nu)_{x}=0,\qquad u_{t}+uu_{x}+n_{x}=0,$ (46)

which represent the dispersionless limit of the NLS equation, can be
transformed to the dispersionless equations (17) by means of substitutions
[23, 14]

$u=2vw,\qquad n=(1-w^{2})(\widetilde{\rho}-v^{2}).$ (47)

Consequently, these substitutions cast the expressions for the canonical
momentum and Hamiltonian of the NLS soliton to the expressions (44) and (45),
respectively. This observation allows us to simplify derivation of the Newton
equation for the magnetic soliton motion.

##  V Newton equation

In Hamiltonian mechanics, we have to solve Eq. (44) with respect to $\dot{x}$
and to substitute this function of $x$ and $p$ into Eq. (45). Then the
Hamilton equations (36) and (39) govern the soliton’s motion. It is impossible
to perform such an elimination of $\dot{x}$ explicitly, therefore it is more
convenient to pass from the Hamilton equations to the Newton-like equation for
the soliton’s acceleration $\ddot{x}$. Actually, this transformation was
already put in practice for the NLS soliton in Ref. [13] and we can use the
formulas obtained there for writing the Newton equation for the magnetic
soliton with help of the substitutions (47). The Newton equation for the NLS
soliton reads [13]

$2\ddot{x}=n_{x}+(u+\dot{x})u_{x}+2u_{t}+\frac{n_{t}+(nu)_{x}}{\sqrt{n-(\dot{x}-u)^{2}}}\arccos\frac{\dot{x}-u}{\sqrt{n}}.$
(48)

Consequently, the substitutions (47) yield the Newton equation for the
magnetic soliton in the following form

$\begin{split}&2\ddot{x}=\left[(1-w^{2})(\widetilde{\rho}-v^{2})\right]_{x}+2(2vw+\dot{x})(vw)_{x}+4(vw)_{t}\\\
&+\frac{\left[(1-w^{2})(\widetilde{\rho}-v^{2})\right]_{t}+\left[2(1-w^{2})(\widetilde{\rho}-v^{2})vw\right]_{x}}{2\sqrt{(1-w^{2})(\widetilde{\rho}-v^{2})-(\dot{x}-2vw)^{2}}}\arccos\frac{\dot{x}-2vw}{\sqrt{(1-w^{2})(\widetilde{\rho}-v^{2})}}.\end{split}$
(49)

As one can see, in this theory the soliton’s dynamic is governed by the
gradients of the background distributions in agreement with our assumption
that the soliton’s width is much smaller than the characteristic length at
which the background changes, so the soliton’s profile is only determined by
the soliton’s velocity and the local values of the background variables. These
variables evolve according to the dispersionless equations obtained from Eqs.
(14) by neglecting the dispersion term in the second equation,

$w_{t}=\left[v(1-w^{2})\right]_{x},\qquad
v_{t}=\left[w(\widetilde{\rho}-v^{2})\right]_{x}+V_{1,x}-V_{2,x}.$ (50)

In derivation of expressions for the canonical momentum and the Hamiltonian,
we assumed that the energy and momentum do not depend on the gradients of the
potentials, so we used Eqs. (17). Now, in discussion of the large scale
dynamics of the background, we have to take into account the gradients of the
potentials $V_{1,x}$ and $V_{2,x}$ as well as slow $x$-dependence of the mean
total density in the variable $\widetilde{\rho}=2\delta\rho_{0}$, where
$\rho_{0}(x)$ is a function of $x$ determined by Eq. (6). Thus, we can
simplify Eq. (49) by exclusion of $w_{t}$ and $v_{t}$ with the use of Eqs.
(50) and after some manipulations we get

$\begin{split}\frac{d^{2}x}{dt^{2}}=\frac{d(vw)}{dt}+\frac{1+w^{2}}{2}\widetilde{\rho}_{x}+w(V_{1}-V_{2})_{x}-\frac{v(1-w^{2})(V_{1}-V_{2})_{x}}{\sqrt{(1-w^{2})(\widetilde{\rho}-v^{2})-(\dot{x}-2vw)^{2}}}\arccos\frac{\dot{x}-2vw}{\sqrt{(1-w^{2})(\widetilde{\rho}-v^{2})}}.\end{split}$
(51)

Let us consider some particular cases.

If $V_{1}=V_{2}=0$ and, consequently, $\widetilde{\rho}=\mathrm{const}$, this
equation reduces to $d^{2}x/dt^{2}=d(vw)/dt$ so its evident integration
reproduces to the relationship (26). This equation is very convenient for
discussions of the soliton’s motion along non-uniform evolving background
waves which are not subject to action of external potentials. For example, if
the soliton moves along a self-similar rarefaction wave evolved from and
initial discontinuity [28, 19], then of the Riemann invariants (20) is
constant. Let it be $r_{-}=r_{-}^{0}=\mathrm{const}$, then one the other
invariant $r_{+}$ is determined by the equation $v_{+}=x/t$, so these two
equations give $r_{-}+v_{+}=3vw=x/t+r_{-}^{0}$. Equation

$\frac{dx}{dt}=\frac{1}{3}\left(\frac{x}{t}+r_{-}^{0}\right)+q=\frac{1}{3}\frac{x}{t}+q_{1},$
(52)

where $q_{1}=q+r_{-}^{0}/3$ is determined by the parameter $r_{-}^{0}$ of the
background flow and by the initial soliton’s velocity at the moment $t=t_{0}$.
This equations is readily integrated to give

$x(t)=\frac{3}{2}q_{1}t+Ct^{1/3},$ (53)

where the integration constant $C$ is determined by the initial soliton’s
coordinate $x_{0}=x(t_{0})$.

Equation (51) greatly simplifies, if there is no relative motion of the
components ($v=0$),

$\frac{d^{2}x}{dt^{2}}=\frac{1+w^{2}}{2}\widetilde{\rho}_{x}+w(V_{1}-V_{2})_{x},$
(54)

where $\widetilde{\rho}(x)=2\delta g\,\rho_{0}(x)=(\delta g/g)(\mu-
V_{1}(x)-V_{2}(x))$ (see Eq. (6)). If the trap potential is the same for both
components, $V_{1}=V_{2}=V_{ext}$, then $w=\mathrm{const}$ and Eq. (54)
reduces to

$\frac{d^{2}x}{dt^{2}}=-\frac{(1+w^{2})\delta g}{g}\frac{dV_{ext}}{dx}.$

Here we should remember that in these calculations we used the space variable
$x^{\prime}=\sqrt{2}x$ (see text below Eq. (10)), so returning to the usual
coordinate $x$, we write this equation as

$\frac{d^{2}x}{dt^{2}}=-\frac{(1+w^{2})\delta g}{2g}\frac{dV_{ext}}{dx}.$ (55)

Multiplication by $\dot{x}$ and integration yield the integral of motion

$\frac{1}{2}\left(\frac{dx}{dt}\right)^{2}+\frac{(1+w^{2})\delta
g}{2g}V_{ext}=\mathrm{const}.$ (56)

This equation determines the soliton’s path along a stationary profile of
condensate’s density $\rho_{0}(x)$ confined in the potential $V_{ext}(x)$. In
a linear potential $V_{ext}=\eta x$ the soliton moves according to the
“Galileo law”

$x=x_{0}-\frac{(1+w^{2})\delta g}{4g}t^{2}.$ (57)

In a harmonic trap $V_{ext}=(\omega_{0}^{2}/2)x^{2}$ the soliton oscillates
with the frequency

$\omega=\sqrt{\frac{(1+w^{2})\delta g}{2g}}\omega_{0}.$ (58)

In case of equal densities of the components ($w=0$), we get

$\omega=\sqrt{\frac{\delta g}{2g}}{\omega_{0}}.$ (59)

It differs from the frequency $\omega=\omega_{0}/\sqrt{2}$ of oscillations of
a dark soliton in a one-component condensate [8, 10] by the factor
$\sqrt{\delta g/g}$.

If $V_{1}=-V_{2}=V_{ext}(x)$ [29], then the total density is constant,
$\rho_{1}+\rho_{2}=\rho_{0}=\mathrm{const}$, as well as
$\widetilde{\rho}=2\delta g\rho_{0}$, but the partial densities are non-
uniform and the second Eq. (50) yields the distribution of $w$,

$w(x)=-\frac{1}{\delta g\,\rho_{0}}V_{ext}(x),$ (60)

where we have chosen the integration constant in such a way that $w=0$ at the
point with $V_{ext}=0$. In this case Eq. (54) reduces with account of the
replacement $x\to\sqrt{2}x$ to

$\frac{d^{2}x}{dt^{2}}=-\frac{1}{2\delta g\rho_{0}}\frac{d(V_{ext})^{2}}{dx},$
(61)

and this equation has an integral of motion

$\frac{1}{2}\left(\frac{dx}{dt}\right)^{2}+\frac{1}{2\delta
g\rho_{0}}V_{ext}^{2}=\mathrm{const}.$ (62)

If $V_{ext}=\eta x$, then soliton oscillates with the frequency

$\omega=\frac{\eta}{\sqrt{\delta g\,\rho_{0}}}.$ (63)

##  VI Conclusion

The method used in this paper is based on the generalization of Stokes’ remark
[12] that soliton’s tails propagate with soliton’s velocity and at the same
time satisfy the linearized equations for small-amplitude waves. Therefore,
the soliton’s velocity is related with the soliton’s inverse half-width by the
dispersion relation continued to the region of complex wave numbers. Validity
of this remark is confirmed by particular soliton solutions of a number of
concrete nonlinear wave equations. It remains also correct for solitons
propagating along non-uniform background, as one can see from the coincidence
of the equation for soliton’s motion with the corresponding results obtained
in the soliton limit of the Whitham modulation equations in completely
integrable cases when the Whitham equations are available in the Riemann
diagonal form. Extension of the obtained by this method equations for solitons
motion to non-integrable situations allows one to reproduce some known results
and to obtain the new ones.

This approach is not limited to the problem of solitons motion. Another its
interesting application is to the quasi-classical method of calculation of
velocities of solitons produced from an intensive initial pulse [30]. In
completely integrable cases, it easily reproduces the known results following
from the WKB method applied to the spectral problem associated with the
concrete equation in framework of the inverse scattering transform method.
What is more important, it can be applied to not completely integrable
equations and gives the results which agree very well with the results of
numerical simulations.

At last, we should mention that our approach allowed us to connect the
condition of preservation of Hamiltonian mechanics of propagation of high-
frequency wave packets along a smooth large-scale dispersionless flow with the
quasi-classical limit of the Lax pair for equations integrable in framework of
AKNS scheme [30, 25]. This observation sheds new light on conditions of
complete integrability of nonlinear wave equations.

One may hope that the suggested method can lead to many other interesting
results.

## Acknowledgments

I am grateful to N. Pavloff and D.V. Shaykin for useful discussions. This
research is funded by the research project FFUU-2024-0003 of the Institute of
Spectroscopy of the Russian Academy of Sciences (Sections 2,3) and by the RSF
grant number 19-72-30028 (Sections 4,5).

## References

  * [1] D. J. Kaup, A Perturbation Expansion for the Zakharov-Shabat Inverse Scattering Transform, SIAM J. Appl. Math. 31, 121-133 (1976).
  * [2] V. I. Karpman, E. M. Maslov, Perturbation theory for solitons, Zh. Eksp. Teor. Fiz. 73, 537 (1977) [Sov. Phys. JETP, 46, 281 (1977)].
  * [3] V. I. Karpman, E. M. Maslov, Structure of shelves created in solitons by perturbations, Zh. Eksp. Teor. Fiz. 75, 504 (1978) [Sov. Phys. JETP, 48, 252 (1978)].
  * [4] Y. Kodama and M.J. Ablowitz, Perturbations of solitons and solitary waves, Stud. Appl. Math., 64, 225-245 (1981).
  * [5] Y. S. Kivshar and B. A. Malomed, Dynamics of solitons in nearly integrable systems, Rev. Mod. Phys., 61, 763 (1989).
  * [6] L. Ostrovsky, K. Gorshkov, Perturbation Theories for Nonlinear Waves, in Nonlinear Science at the Dawn of the 21st Century, Eds. P. L. Christiansen, M. P. Sørensen, A. C. Scott, p. 47, (Berlin, Springer, 2000).
  * [7] Y. S. Kivshar, X. P. Yang, Perturbation-induced dynamics of dark solitons, Phys. Rev. E 49, 1657-1670 (1994).
  * [8] Th. Busch and J. R. Anglin, Motion of Dark Solitons in Trapped Bose-Einstein Condensates, Phys. Rev. Lett. 84, 2298 (2000).
  * [9] M. J. Ablowitz, S. D. Nixon, T. P. Horikis, D. J. Frantzeskakis, Perturbations of dark soliton, Proc. Roy. Soc. London, A 467, 2597-2621 (2011).
  * [10] V. V. Konotop and L. P. Pitaevskii, Landau Dynamics of a Grey Soliton in a Trapped Condensate, Phys. Rev. Lett. 93, 240403 (2004).
  * [11] A. M. Kamchatnov and D. V. Shaykin, Propagation of generalized Korteweg-de Vries solitons along large scale waves, Phys. Rev. E 108, 054205 (2023).
  * [12] G. G. Stokes, Mathematical and Physical Papers, Vol. V, p. 163 (Cambridge University Press, Cambridge, 1905).
  * [13] S. K. Ivanov and A. M. Kamchatnov, Motion of dark solitons in a non-uniform flow of Bose-Einstein condensate, Chaos, 32, 113142 (2022).
  * [14] A. M. Kamchatnov, Hamilton theory of NLS equation soliton motion, Theor. Math. Phys. (to be published); arXiv:2310.04001.
  * [15] A. M. Kamchatnov and D. V. Shaykin, Propagation of dark solitons of DNLS equation along a large-scale background, arXiv:2402.12776.
  * [16] C. Qu, L. P. Pitaevskii, and S. Stringari, Magnetic Solitons in a Binary Bose-Einstein Condensate, Phys. Rev. Lett. 116, 160402 (2016).
  * [17] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation and Superfluidity, (Oxford, UP, 2016)
  * [18] E. J. Mueller, Spin textures in slowly rotating Bose-Einstein condensates, Phys. Rev. A 69, 033606 (2004).
  * [19] T. Congy, A. M. Kamchatnov, N. Pavloff, Dispersive hydrodynamics of nonlinear polarization waves in two-component Bose-Einstein condensates, SciPost Phys. 1, 006 (2016).
  * [20] D. T. Son, M. A. Stephanov, Domain walls of relative phase in two-component Bose-Einstein condensates, Phys. Rev. A 65, 063621 (2002).
  * [21] A. M. Kosevich, B.A. Ivanov, A.S. Kovalev, Magnetic solitons, Phys. Rep. 194, 117-238 (1990).
  * [22] E. Iacocca, M. A. Hoefer, Vortex-antivortex proliferation from an obstacle in thin film ferromagnets, Phys. Rev. B 95, 134409 (2017).
  * [23] L. V. Ovsyannikov, Two-layer “shallow water” model, J. Appl. Mech. Tech. Phys. 20, 127-135 (1979).
  * [24] S. K. Ivanov, A. M. Kamchatnov, T. Congy, and N. Pavloff, Solution of the Riemann problem for polarization waves in a two-component Bose-Einstein condensate, Phys. Rev. E 96, 062202 (2017).
  * [25] A. M. Kamchatnov and D. V. Shaykin, Quasiclassical integrability condition in AKNS scheme, Physica D, 460, 134085 (2024).
  * [26] P. Sprenger, M. A. Hoefer, and G. A. El, Hydrodynamic optical soliton tunneling, Phys. Rev. E 97, 032218 (2018).
  * [27] K. van der Sande, G. A. El, and M. A. Hoefer, Dynamic soliton-mean flow interaction with non-convex flux, J. Fluid Mech. 928, A21 (2021).
  * [28] C. Hamner, J.J. Chang, P. Engels, M.A. Hoefer, Generation of Dark-Bright Soliton Trains in Superfluid-Superfluid Counterflow, Phys. Rev. Lett., 106, 065302 (2011).
  * [29] S. Bresolin, A. Roy, G. Ferrari, A. Recati, and N. Pavloff, Oscillating Solitons and ac Josephson Effect in Ferromagnetic Bose-Bose Mixtures, Phys. Rev. Lett. 130, 220403 (2023).
  * [30] A. M. Kamchatnov, Asymptotic theory of not completely integrable soliton equations, Chaos 33, 093105 (2023).