Hicham Saylani 11institutetext: Laboratoire de matériaux, signaux et systèmes
et modélisation physique,
Faculté des Sciences, Université Ibn Zohr, BP 8106, Cité Dakhla, Agadir,
Morocco 11email<EMAIL_ADDRESS>11email<EMAIL_ADDRESS>
# Identification of melanoma diseases from multispectral dermatological images
using a novel BSS approach

Mustapha Zokay  Hicham Saylani

###### Abstract

In this paper111The paper has been accepted at the MICAD 2022 conference and
is currently being published as a chapter in the Springer book :  Lecture
Notes in Electrical Engineering. we propose a new approach to identify
melanoma diseases by identifying the distribution of its main skin
chromophores (melanin, oxyhemoglobin and deoxyhemoglobin) from multispectral
dermatological images. Based on Blind Source Separation (BSS), our approach
takes into account the shading present in most of the images. Assuming that
the multispectral images have at least 4 spectral bands, it allows to estimate
the distribution of each chromophore in addition to the shading without any a
priori information, contrary to all existing methods that use 3 bands, i.e.
RGB images. Indeed, the fact of neglecting the shading degrades their
performance. To validate our method, we used a database of real multispectral
dermatological images of skin affected by melanoma cancer. To measure our
performance, in addition to the classical criterion of visually analyzing the
estimated distributions with referring to the physiological knowledge of the
disease, we proposed a new criterion that is based on our independence
hypothesis. Using these two criteria, we could see that our approach is very
efficient for the identification of melanoma.

###### keywords:

Multispectral dermatological images, Chromophores, Melanin, Hemoglobin,
Oxyhemoglobin, Deoxyhemoglobin, Shading, Blind Source Separation (BSS),
Melanoma.

##  1 Introduction

The skin is the largest organ in the human body. It contains three main
chromophores which are melanin, oxyhemoglobin and deoxyhemoglobin. The
distribution of these chromophores is of great importance for dermatologists
who use it for the identification and early detection of skin diseases. One of
the techniques that is increasingly used and has proven its effectiveness in
identifying the distribution of chromophores is multispectral imaging [1].
However, the direct use of multispectral images tends to give erroneous
information on the distribution of chromophores. Indeed, the light intensity
reflected by the skin does not only depend on these three chromophores but
also on the geometry of the skin surface, called shading. The multispectral
images obtained at different wavelengths can therefore be considered as
mixtures of four constituents which are the three chromophores and the
shading. Thus, if we consider these constituents as sources, this problem can
be seen as a source separation problem, known as an inverse problem that
belongs to the field of signal processing. The objective of source separation
is to estimate the sources from the sole knowledge of their mixtures. As this
estimation is done without any a priori information neither on the sources nor
on the mixing coefficients, it is called Blind Source Separation (BSS). It is
easy to see that the BSS problem is an ill-posed inverse problem which admits
an infinity of solutions so it is essential to add hypotheses on the sources
and/or on the mixing coefficients, which gave rise to 3 main families of BSS
methods which are the Independent Component Analysis (ICA), the Sparse
Component Analysis (SCA) and the Non-negative Matrix Factorization (NMF) (see
[2] for more details). During the last decade, the use of BSS methods for non-
invasive identification of chromophore distributions has attracted the
interest of several researchers [3, 4, 5, 6, 7, 8, 9] who all adopted the same
mixing model. Knowing that most of these researchers were interested in RGB
images and thus in 3 wavelengths $\lambda_{i}$, $i\in${1,2,3} which represent
respectively the central wavelengths of the Blue, Green and Red bands, and if
we note $j\in${$1,2,3$} the index related to the chromophores and which
represents respectively melanin, oxyhemoglobin and deoxyhemoglobin, then the
classical mixing model is written:

$I_{\lambda_{i}}(\textit{{u}})=\sum_{j=1}^{j=3}a_{ij}.S_{j}(\textit{{u}})+p_{d}(\textit{{u}})+n_{i},\hskip
11.38092pti\in\\{1,2,3\\},\vspace*{-0.5cm}$ (1)

where:

  * $\bullet$

$I_{\lambda_{i}}(\textit{{u}})$ is the logarithm inverse of the reflectance
detected by the camera, at wavelength $\lambda_{i}$, at the pixel of
coordinates $(x,y)=\textit{{u}}$,

  * $\bullet$

$S_{j}(\textit{{u}})$ represents the chromophore of index $j$,

  * $\bullet$

$a_{ij}$ represents the mixing coefficient which depends on the molar
absorption coefficient of the chromophore $S_{j}(\textit{{u}})$ and the light
penetration depth in the skin at wavelength $\lambda_{i}$,

  * $\bullet$

$p_{d}(\textit{{u}})$ represents the shading variation in the image,

  * $\bullet$

$n_{i}$ represents the characteristics of the sensor.

This mixing model has been adopted by all existing BSS methods [3, 4, 5, 6, 7,
8], but they differ in the assumptions made about the chromophores and the
procedure followed to reach the final objective, so that these methods can be
classified into two main families. The first family includes the BSS methods
that consider oxyhemoglobin and deoxyhemoglobin as a single source called
hemoglobin [3, 4, 5, 10, 9]. Indeed, in [5, 10], the authors applied Principal
Component Analysis (PCA) and then ICA on an RGB image of the skin, assuming
that the shading is constant throughout the image and that melanin and
hemoglobin are independent. In [11], Madooei et al. proposed a new 2-D color
chromaticity to eliminate shading using the geometric mean color, and then
estimated the distributions of the two chomophores using ICA. In [9], Gong et
al. proposed to estimate the distributions of both chromophores from an RGB
image using NMF. In [3], Galeano et al. proposed to use a neural network-based
system and then applied NMF to separate the melanin from the hemoglobin, but
in their model they neglected the shading.
The second family includes the BSS methods which are based on a priori
information on the absorption spectra of the three chromophores and the light
penetration depth in the skin [6, 7]. the authors proposed to remove the
shading and the specular reflection from the RGB image using respectively
white paper and polarizers, then they relied on knowledge of absorption
coefficients and estimates the light penetration depth in order to deduce an
empirical mixing matrix to extract the distributions of the three
chromophores. It should be noted, that the weakness of this family lies in the
level of accurate estimation of light penetration depth into the skin.

In this paper, we propose a new method based on BSS to estimate the
distribution of all the three main chromophores separately, in addition to the
shading distribution which we consider as a full-fledged source, unlike all
existing methods. Based on more realistic assumptions and applying to
multispectral images with at least 4 spectral bands, our new method exploits
the intrinsic properties of each chromophore. To validate our method we use a
database of real multispectral dermatological images of skin affected by
melanoma cancer disease [12]. For the performance measurement, in addition to
the classical criterion which is based on the visual analysis of the estimated
distributions of the three chromophores, we propose in this paper a new
numerical criterion which is based on the measure of independence between the
estimated distributions of melanin and hemoglobin. The rest of this paper is
organized as follows. Section 2 presents our new method for estimating the
distribution of the three chromophores and shading. Section 3 presents the
results of the tests performed followed by a last section dedicated to a
conclusion and perspectives for our work.

##  2 Proposed method

The first key idea of our new method is to consider shading as a full-fledged
source, in addition to the three sources of interest, which allows us to avoid
the unrealistic assumption made by most existing methods that its contribution
is the same at all wavelengths222Indeed, we found from the experimental curve
obtained by PCA used in [13] that the contribution of shading ($p_{d}$) is not
equal to $1$ in all mixtures.. However, as the number of sources involved
becomes equal to 4 we are interested in this paper for the case where we have
multispectral images with at least 4 bands, said case determined333This is the
case where we have as many mixtures as sources. On the other hand, in the case
where we have more mixtures than sources, called the over-determined case, we
can easily return to the determined case by applying a PCA.. From equation
(1), we can see that the term $n_{i}$ does not give any information on
chromophore distributions. Thus, as in [5], we begin by eliminating it from
our mixtures based on the following hypothesis (H1).

(H1) :

There is at least one pixel in the 4 images where the concentrations of the
three chromophores and shading are all zero, i.e.:

$\exists\,\,\textit{{u}}\,\,/\,\,n_{i}=min(I_{\lambda_{i}}(\textit{{u}})).$
(2)

So, the new mixture model is written:

$X_{i}(\textit{{u}})=\sum_{j=1}^{j=3}a_{ij}\cdot
S_{j}(\textit{{u}})+a_{i4}\cdot S_{4}(\textit{{u}}),\hskip
11.38092pti\in[1,4]\vspace*{-0.3cm}$ (3)

where $\,\,X_{i}(\textit{{u}})=I_{\lambda_{i}}(\textit{{u}})-n_{i}$,
$S_{4}(\textit{{u}})=p_{d}(\textit{{u}})$ and $a_{i4}\in\mathbb{R}$.

As for all BSS methods, we generate new 1D mixtures (vectors), which we note
$X_{i}(v)$, from the 2D mixtures (images) $X_{i}(\textit{{u}})$ by
concatenating the rows of the latter. We then have:

$X_{i}(v)=\text{vec}(X_{i}(\textit{{u}}))=\sum_{j=1}^{j=4}a_{ij}\cdot
S_{j}(v),\hskip 11.38092pti\in[1,4]$ (4)

The second key idea of our method is to treat mixtures in two steps, unlike
all existing methods that treat all mixtures at the same time. Indeed, we
start by treating only two mixtures that contain only melanin and shading in
order to separate them first, and then we eliminate their contribution from
the other two mixtures remaining to keep only oxyhemoglobin and
deoxyhemoglobin. These last two chromophores are then separated in a last
step. These three steps of our method are detailed below.

### Step 1: Separation of sources $S_{1}(v)$ and $S_{4}(v)$

In this step we exploit the properties of each chromophore concerning spectral
absorption as a function of the wavelength. Indeed, based on data published in
[14], we found that light absorption at wavelengths greater than $620\,nm$ is
dominated by melanin, so the absorption coefficients $a_{ij}$ of oxyhemoglobin
and deoxyhemoglobin are all negligible, i.e. we have $a_{32}=a_{33}=0$ and
$a_{42}=a_{43}=0$. Thus, the mixtures corresponding to the red and infrared
bands can be re-written as follows:

$\left\\{\begin{array}[]{ll}X_{3}(v)=a_{31}\cdot S_{1}(v)+a_{34}\cdot
S_{4}(v)\\\ X_{4}(v)=a_{41}\cdot S_{1}(v)+a_{44}\cdot
S_{4}(v)\end{array}\right.$ (5)

We can re-write the equation system (5) in a matrix form as follows:

${\bf X}(v)={\bf A}\cdot{\bf S}(v),$ (6)

where ${\bf X}(v)=[X_{3}(v),X_{4}(v)]^{T}$, ${\bf
S}(v)=[S_{1}(v),S_{4}(v)]^{T}$ and ${\bf A}=\begin{pmatrix}a_{31}&a_{34}\\\
a_{41}&a_{44}\\\ \end{pmatrix}$.

In this first step we assumed that the two sources $S_{1}(v)$ and $S_{4}(v)$
are independent, in which case we can use one of the ICA methods to separate
them. We have opted here for the AMUSE method [15] for its simplicity since it
exploits only the second order statistics of the signals. Indeed, the working
hypotheses of this method are the following.

(H2) :

The sources $S_{j}(v)$ are auto-correlated and mutually uncorrelated, i.e. :

$\forall\,\,\tau,\,\hskip 8.5359pt\left\\{\begin{array}[]{ll}\hskip
2.84544ptE[S_{j}(v)\cdot S_{j}(v-\tau)]\neq 0,\,\,\,\,j\in\\{1,4\\}\\\ \hskip
2.84544ptE[S_{1}(v)\cdot S_{4}(v-\tau)]=E[S_{1}(v)]\cdot
E[S_{4}(v-\tau)]\end{array}\right.$ (7)

(H3) :

The condition of identifiability for the method is verified, i.e. :

$\exists\,\,\tau\neq 0\,\,/\hskip 8.5359pt\frac{E[S_{1}(v)\cdot
S_{1}(v-\tau)]}{E[S_{1}^{2}(v)]}\neq\frac{E[S_{4}(v)\cdot
S_{4}(v-\tau)]}{E[S_{4}^{2}(v)]}.$ (8)

The method AMUSE allows us to estimate the separation matrix ${\bf A}^{-1}$ to
a permutation matrix ${\bf P}$ and a diagonal matrix ${\bf D}$ [15]. By noting
this matrix ${\bf C}={\bf P}{\bf D}{\bf A}^{-1}$, we obtain finally the source
matrix ${\bf S}(v)$ with the same indeterminations as follows:

${\bf C}\cdot{\bf X}(v)=\left({\bf PDA}^{-1}\right)\cdot\left({\bf
AS}(v)\right)={\bf PD}{\bf S}(v).$ (9)

Indeed, by noting ${\bf PD}{\bf S}(v)={\bf Y}(v)=[Y_{1}(v),Y_{4}(v)]$ and
omitting the permutation444Indeed, the permutation matrix $\bf P$ can be
identified based on the visual analysis since the shading source $S_{4}(v)$ is
easily differentiable. we have:

$Y_{j}(v)=\alpha_{j}\cdot S_{j}(v),\hskip 5.69046ptj=1,4,$ (10)

where the $\alpha_{j}$ are the elements constituting the diagonal of the
matrix ${\bf D}$.

### Step 2: Removal of $S_{1}(v)$ and $S_{4}(v)$ sources from mixtures

The goal of this step is to eliminate the contributions of the sources
estimated $S_{1}(v)$ and $S_{4}(v)$ from the mixtures $X_{1}(v)$ and
$X_{2}(v)$. For this we exploit the following independence hypothesis.

(H4) :

$S_{i}(v)$ and $S_{j}(v)$ are mutually uncorrelated instantaneously, for
$i\in\left\\{2,3\right\\}$ and $j\in\left\\{1,4\right\\}$, i.e.:

$E[\tilde{S}_{i}(v)\cdot\tilde{S}_{j}(v)]=E[\tilde{S}_{i}(v)]\cdot
E[\tilde{S}_{j}(v)]=0,\hskip
2.84544pt\forall(i,j)\in\left\\{2,3\right\\}\times\left\\{1,4\right\\}$ (11)

where $\tilde{S}_{i}(v)$ are the centered versions555i.e. :
$\tilde{S}_{i}(v)=S_{i}(v)-E[S_{i}(v)]$. of the sources $S_{i}(v)$. By
denoting respectively $\tilde{X_{i}}(v)$ and $\tilde{Y_{j}}(v)$ the centered
versions of the signals $X_{i}(v)$ and $Y_{j}(v)$, and by exploiting the
relations (10) and (11) we can write:

$\displaystyle E[\tilde{X_{i}}(v)\cdot\tilde{Y}_{j}(v)]$ $\displaystyle=$
$\displaystyle
E\left[\left(\sum_{k=1}^{k=4}a_{ik}\cdot\tilde{S}_{k}(v)\right)\cdot\left(\alpha_{j}\cdot\tilde{S}_{j}(v)\right)\right]$
(12) $\displaystyle=$ $\displaystyle a_{ij}\cdot\alpha_{j}\cdot
E[\tilde{S}_{j}^{2}(v)]$ (13)

On the other hand, according to the relation (10) we have:

$E[\tilde{Y}_{j}^{2}(v)]=\alpha_{j}^{2}\cdot E[\tilde{S_{j}^{2}}(v)],\hskip
11.38092ptj=1,4.$ (14)

Thus, by exploiting the relations (13) and (14) we can generate two new
mixtures $Z_{i}(v)(i=1,2)$ which contain only the sources $S_{2}(v)$ and
$S_{3}(v)$ as follows:

$\displaystyle Z_{i}(v)$ $\displaystyle=$ $\displaystyle
X_{i}(v)-\frac{E[\tilde{X_{i}}(v)\cdot\tilde{Y_{1}}(v)]}{E[\tilde{Y_{1}^{2}}(v)]}\cdot
Y_{1}(v)-\frac{E[\tilde{X_{i}}(v)\cdot\tilde{Y_{4}}(v)]}{E[\tilde{Y_{4}^{2}}(v)]}\cdot
Y_{4}(v)$ (15) $\displaystyle=$ $\displaystyle
X_{i}(v)-\frac{a_{i1}}{\alpha_{1}}\cdot(\alpha_{1}.\cdot
S_{1}(v))-\frac{a_{i4}}{\alpha_{4}}\cdot(\alpha_{4}\cdot S_{4}(v))$ (16)
$\displaystyle=$ $\displaystyle a_{i2}\cdot S_{2}(v)+a_{i3}\cdot S_{3}(v)$
(17)

### Step 3 : Separation of Sources $S_{2}(v)$ and $S_{3}(v)$

The goal of this step is to separate the remaining $S_{2}(v)$ and $S_{3}(v)$
sources that represent oxyhemoglobin and deoxyhemoglobin respectively, and
this time by treating the new mixtures $Z_{1}(v)$ and $Z_{2}(v)$ provided by
the previous step. As these two sources are dependent, which is why most of
the existing methods fail to separate them, we have opted in this paper for a
new solution based on the exploitation of their sparsity666A source is said to
be sparse in a given representation domain if some of its samples are zero in
this domain.. Here is our working hypothesis for this step.

(H5) :

For each source, there is at least one spatial area over which it is active
while the other source is inactive, i.e.:

$\left\\{\begin{array}[]{ll}\exists\,V_{1}\,/\,\forall\,v\in
V_{1},\,\,{S_{2}(v)=0}\,\,\,\,\text{and}\,\,\,\,{S_{3}(v)\neq 0}\\\
\exists\,V_{2}\,/\,\forall\,v\in
V_{2},\,\,{S_{3}(v)=0}\,\,\,\,\text{and}\,\,\,\,{S_{2}(v)\neq
0}\end{array}\right.$ (18)

There are several BSS methods that exploit this sparsity assumption to achieve
separation. We opted here for the TEMPROM method proposed in[16] for its
simplicity. This method consists of identifying the single-source areas
$V_{1}$ and $V_{2}$ at first, and then calculating the ratio between the
mixtures $Z_{2}(v)$ and $Z_{1}(v)$ on these areas at a second time, which
ultimately allows to estimate the matrix of mixture involved. Indeed, by
exploiting the two equations of the system (18) and the relation (17) we
obtain:

$\left\\{\begin{array}[]{ll}\forall\,v\in
V_{1},\,\,\,\frac{Z_{2}(v)}{Z_{1}(v)}=\frac{a_{23}\cdot S_{3}(v)}{a_{13}\cdot
S_{3}(v)}=\frac{a_{23}}{a_{13}}=r_{1}\\\ \forall\,v\in
V_{2},\,\,\,\frac{Z_{2}(v)}{Z_{1}(v)}=\frac{a_{22}\cdot S_{2}(v)}{a_{12}\cdot
S_{2}(v)}=\frac{a_{22}}{a_{12}}=r_{2}\end{array}\right.$ (19)

Finally, by exploiting the relations (17) and (19) we obtain:

$\left\\{\begin{array}[]{ll}r_{1}\cdot Z_{1}(v)-Z_{2}(v)=(r_{1}\cdot
a_{12}-a_{22})\cdot S_{2}(v)=\alpha_{2}\cdot S_{2}(v)=Y_{2}(v)\\\ r_{2}\cdot
Z_{1}(v)-Z_{2}(v)=(r_{2}\cdot a_{13}-a_{23})\cdot S_{3}(v)=\alpha_{3}\cdot
S_{3}(v)=Y_{3}(v)\end{array}\right.$ (20)

where $\alpha_{2}$ and $\alpha_{3}$ are scalars.

##  3 Results

In this section, we measure the performances of our method using a database of
melanoma skin cancer patients which is an open access database [12]. It is
recalled that our method makes it possible to estimate the distributions of
oxyhemoglobin and deoxyhemoglobin separately in addition to those of melanin
and shading contrary to all existing methods (as mentioned in the
introduction). So, in absolute terms we cannot compare our performance with
any of these methods. However, since most of these methods estimate melanin
and hemoglobin (which is a mixture of oxyhemoglobin and deoxyhemoglobin), we
can limit ourselves in the comparison to these two chromophores. For this, we
opted for a comparison with the methods proposed in [9, 11] because they are
accessible for testing, unlike most of the existing methods [4, 6, 7]. As for
the performance measurement criteria, in addition to the classical visual
criterion which is a subjective criterion, we propose in this paper a new
numerical criterion. In fact, for dermatologists, melanoma is characterized by
a high distribution of melanin, an average distribution of deoxyhemoglobin and
a very low distribution of oxyhemoglobin compared to healthy skin [17]. To
support this subjective criterion, which is the most used in the literature
[9, 4, 6, 7], we use a second numerical criterion, which allows us to check to
what extent the independence hypothesis we have assumed is satisfied by the
estimated chromophores, since this hypothesis is the most used by researchers
[11, 5, 10]. We then define our new performance measurement criterion, which
we note as $C_{Ind}$, by exploiting the assumption of statistical independence
at order 4 between melanin and hemoglobin, as follows:

$C_{Ind}=\frac{1}{2}\left(C_{13}+C_{31}\right),$ (21)

where:

$C_{13}=20\,log_{10}\left(\left|E[\tilde{Y}_{1}(u)\tilde{Z}_{1}(u)^{3}]\right|^{-1}\right)\hskip
5.69046pt\text{and}\hskip
5.69046ptC_{31}=20\,log_{10}\left(\left|E[\tilde{Y}_{1}(u)^{3}\tilde{Z}_{1}(u)]\right|^{-1}\right)$
(22)

We recall that the estimates $\tilde{Y}_{1}(u)$ and $\tilde{Z}_{1}(u)$ are
respectively the centered versions of the melanin and hemoglobin
distributions.

The database contains 30 multispectral images [12]. Knowing that we have
tested our method on some images and that the results obtained are similar, we
limit ourselves here to present the results of a single melanoma image, due to
lack of space. This image is shown in Figure 1.

Figure 1: Treated multispectral dermatological image.

For the methods [11, 9], we used respectively the algorithm FastICA [18] and
MU [19]. Two-dimensional (2D) and three-dimensional (3D) representations of
the estimated distributions of each of the three chromophores in addition to
the shading using our method are given in Figure 2. These 2D and 3D
representations are grouped by column for each chromophore. The 2D and 3D
representations of the estimated melanin and hemoglobin distributions using
the methods [9, 11] are grouped in Figure 3.

|  |  |
---|---|---|---
|  |  |
| Oxyhemoglobin |  | Shading
Melanin |  | Deoxyhemoglobin |
Figure 2: Estimated distributions (2D and 3D) of the three chromophores and
shading using our method.

|  |  |
---|---|---|---
|  |  |
| Hemoglobin |  | Hemoglobin
Melanin |  | Melanin |

(a) Method [9] |  |  |  (b) Method [11]
---|---|---|---

Figure 3: Estimated distributions (2D and 3D) of melanin and hemoglobin using:
(a) Method [9], (b) Method [11].

Figure 2 shows the good performance of our method for the estimation of the
distribution of each of the three chromophores in addition to the shading.
Indeed, we see that in the melanoma area we have a high distribution of
melanin and a relatively high distribution of deoxyhemoglobin compared to
oxyhemoglobin as shown in the 3D representations of these chromophores. This
is fully consistent with the physiological knowledge that characterizes
melanoma, since naevus disease, which generally has a melanin distribution
similar to that of melanoma, is instead characterized by an oxyhemoglobin
distribution that is close to that of deoxyhemoglobin [17]. On the other hand,
from Figure 3, we see that the melanin distribution estimated by the methods
[9, 11] is similar to that estimated by our method, which means that these two
methods [9, 11] also suspects melanoma. Nevertheless, its estimated hemoglobin
distribution (which can be seen as a mixture of oxyhemoglobin,
deoxyhemoglobin, and shading) does not allow to decide on the nature of the
disease, since the distinction between melanoma and nevus can be made only by
estimating the distributions of oxyhemoglobin and deoxyhemoglobin separately.
All these findings, which were made on the basis of physiological knowledge of
the melanoma disease, are in perfect agreement with the results obtained using
our new numerical criterion, denoted $C_{Ind}$ and defined by the relation
(21), which are presented in the table 1. In this table we provide the mean
and standard deviation of $C_{Ind}$, denoted $\overline{C}_{Ind}$ and
$\sigma$, obtained on 10 different images from the database [12].

| Method [9] | Method [11] | Our method
---|---|---|---
$\overline{C}_{Ind}\,(dB)$ | -2.05 | 27.45 | 38.44
$\sigma\,(dB)$ | 3.66 | 6.20 | 4.33
Table 1: Mean and standard deviation of $C_{Ind}$ in $dB$.

From the table 1, we see that our method is much better than the other methods
since the value obtained for $\overline{C}_{Ind}$ using our method is larger
than that obtained using the methods [9, 11]. We also find that the [9] method
has poor performance compared to our method and the method [11] and this can
be explained by the infinite solution problem posed by NMF. Since this
criterion is based on a measure of independence, we deduce that our method
provides estimates of melanin and hemoglobin distributions in output that are
much more independent than those provided by the methods [9, 11].

##  4 Conclusion

In this paper we have proposed a new approach which aims to identify the
melanoma diseases by identifying the distribution of its main skin
chromophores (melanin, oxyhemoglobin and deoxyhemoglobin) from multispectral
dermatological images. The key idea of our approach is to take into account
the shading, considering it as a full-fledged source, and the three
chromophores, which leads to an instantaneous linear mixture model with four
sources rather than two or three sources, as is the case for all existing
methods. The results of all the tests carried out, using a database of real
multispectral images of skin affected by a skin cancer of the type melanoma,
have shown that our approach is very efficient using the classical criterion
based on visual analysis than our new independence criterion. In terms of
perspective, it would be interesting to test our method on other multispectral
dermatological image databases of other skin diseases.

## References

  * [1] Jacques, S.L., Samatham, R., Choudhury, N.: Rapid spectral analysis for spectral imaging. Biomedical optics express 1(1), 157–164 (2010)
  * [2] Comon, P., Jutten, C.: Handbook of Blind Source Separation, Independent Component Analysis and Applications (02 2010)
  * [3] Jolivot, R., Marzani, F., et al.: Quantification of melanin and hemoglobin in humain skin from multispectral image acquisition: use of a neuronal network combined to a non-negative matrix factorization. Applied and Computational Mathematics, special issue on Applied Artificial Intelligence and Soft Computing 11(2), 257–270 (2012)
  * [4] Liu, Z., Zerubia, J.: Melanin and hemoglobin identification for skin disease analysis. In: 2013 2nd IAPR Asian Conference on Pattern Recognition. pp. 145–149. IEEE (2013)
  * [5] Ojima, N., Akazaki, S., Hori, K., Tsumura, N., Miyake, Y.: Application of image-based skin chromophore analysis to cosmetics. Journal of Imaging Science and Technology 48(3), 222–226 (2004)
  * [6] Spigulis, J., Oshina, I.: Snapshot rgb mapping of skin melanin and hemoglobin. Journal of biomedical optics 20(5), 050503 (2015)
  * [7] Spigulis, J., Oshina, I., Berzina, A., Bykov, A.: Smartphone snapshot mapping of skin chromophores under triple-wavelength laser illumination. Journal of Biomedical Optics 22(9), 091508 (2017)
  * [8] Kuzmina, I., Diebele, I., Asare, L., Kempele, A., Abelite, A., Jakovels, D., Spigulis, J.: Multispectral imaging of pigmented and vascular cutaneous malformations: the influence of laser treatment. In: Laser Applications in Life Sciences. vol. 7376, p. 73760J. International Society for Optics and Photonics (2010)
  * [9] Gong, H., Desvignes, M.: Hemoglobin and melanin quantification on skin images. In: International Conference Image Analysis and Recognition. pp. 198–205. Springer (2012)
  * [10] Mitra, J., Jolivot, R., Vabres, P., Marzani, F.S.: Source separation on hyperspectral cube applied to dermatology. In: Medical Imaging 2010: Computer-Aided Diagnosis. vol. 7624, p. 76243I. International Society for Optics and Photonics (2010)
  * [11] Madooei, A., Drew, M.: A Bioinspired Color Representation for Dermoscopy Image Analysis, pp. 23–66 (09 2015)
  * [12] Lézoray, O., Revenu, M., Desvignes, M.: Graph-based skin lesion segmentation of multispectral dermoscopic images. In: International Conference on Image Processing (IEEE). pp. 897–901 (2014)
  * [13] Eguizabal, A., Laughney, A.M., García-Allende, P.B., Krishnaswamy, V., Wells, W.A., Paulsen, K.D., Pogue, B.W., Lopez-Higuera, J.M., Conde, O.M.: Direct identification of breast cancer pathologies using blind separation of label-free localized reflectance measurements. Biomedical optics express 4(7), 1104–1118 (2013)
  * [14] Van Gemert, M., Jacques, S.L., Sterenborg, H., Star, W.: Skin optics. IEEE Transactions on biomedical engineering 36(12), 1146–1154 (1989)
  * [15] Tong, L., Liu, R.W., Soon, V.C., Huang, Y.F.: Indeterminacy and identifiability of blind identification. IEEE Transactions on circuits and systems 38(5), 499–509 (1991)
  * [16] Abrard, F., Deville, Y.: A time–frequency blind signal separation method applicable to underdetermined mixtures of dependent sources. Signal processing 85(7), 1389–1403 (2005)
  * [17] Lihacova, I.: Evaluation of skin oncologic pathologies by multispectral imaging methods. Ph.D. thesis (07 2015)
  * [18] Hyvärinen, A., Oja, E.: A fast fixed-point algorithm for independent component analysis. Neural computation 9(7), 1483–1492 (1997)
  * [19] Lee, D., Seung, H.S.: Algorithms for non-negative matrix factorization. Advances in neural information processing systems 13 (2000)


# Incommensurable matter-wave jets in quasi-1D geometry

Tadej Mežnaršič<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova 39,
SI-1000 Ljubljana, Slovenia  Faculty of mathematics and physics, University of
Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia  Rok Žitko  Jožef Stefan
Institute, Jamova 39, SI-1000 Ljubljana, Slovenia  Faculty of mathematics and
physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia
Katja Gosar  Jožef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia
Faculty of mathematics and physics, University of Ljubljana, Jadranska 19,
SI-1000 Ljubljana, Slovenia  Katja Arh  Jožef Stefan Institute, Jamova 39,
SI-1000 Ljubljana, Slovenia  Faculty of mathematics and physics, University of
Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia  Matevž Jug  Jožef Stefan
Institute, Jamova 39, SI-1000 Ljubljana, Slovenia  Faculty of mathematics and
physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia
Erik Zupanič  Jožef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia
Faculty of natural sciences and engineering, University of Ljubljana,
Aškerčeva 12, SI-1000 Ljubljana, Slovenia  Peter Jeglič<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia  Faculty of
mathematics and physics, University of Ljubljana, Jadranska 19, SI-1000
Ljubljana, Slovenia

###### Abstract

We experimentally show the formation of incommensurable “golden”
$\frac{1+\sqrt{5}}{2}$ matter-wave jets in a Bose-Einstein condensate (BEC)
subjected to single frequency interaction modulation. We study the formation
of higher order jets and the corresponding incommensurable density waves in
quasi one dimensional (1D) geometry with the help of numerical 1D Gross-
Pitaevskii equation simulation. We explore the process of jet formation
experimentally and theoretically for a large range of modulation amplitudes
and frequencies and present a phase diagram for jet formation. In simulation,
for large modulation amplitudes, the distribution of jet velocities forms a
supercontinuum, which emerges from the observed incommensurable matter-wave
jets.

Incommensurability is encountered in aperiodic crystals with incommensurable
phases [1, 2, 3], charge and spin structures [4, 5, 6, 7, 8], twisted moiré
bilayer experiments [9] with Hofstadter butterfly energy structure [10], and
surface layers [11]. In the context of cold atoms, incommensurable optical
lattices have been used to study wavefunction localization phenomena in
disordered potentials within the Aubry-André model [12, 13, 14, 15, 16, 17,
18, 19, 20].

In this Letter we experimentally show the formation of incommensurable density
waves inside a cesium BEC confined in a quasi one-dimensional (1D) potential
evidenced via the emission of matter-wave jets. Previous studies of matter-
wave jets in two dimensions have revealed the formation of density waves
inside the condensate prior to the emission of jets and explained the process
of formation of higher-order jets [21, 22]. We expand on our previous
experiments with matter-wave jets in a quasi-1D potential [23] by
experimentally and numerically exploring the formation of incommensurable
“golden” jets and other higher-order jets in a quasi-1D BEC under a single-
frequency interaction modulation.

Figure 1:  (a) Absorption images of three different realizations of golden
jets (for modulation frequency 4 kHz, amplitude
$a_{\mathrm{ac}}=90a_{\mathrm{0}}$, $t_{\mathrm{p}}=10\,\mathrm{ms}$, and
$t_{\mathrm{e}}=40\,\mathrm{ms}$). We see the original BEC (J0), first order
jets (J1), second order jets (J2) and two pairs of “golden” jets (G2) (but not
$\sqrt{2}$ D2 jets). Each image shows an area of 1037 µm $\times$ 107 µm. (b)
Schematic of possible jets for increasing number of subsequent collisions
(order, encoded by color).

We start by preparing an almost pure BEC of approximately 4000 cesium atoms in
a crossed dipole trap, from which we release it into a channel with radial
frequency $\omega_{\mathrm{r}}=2\pi\cdot 90\,\mathrm{Hz}$. Simultaneously, we
change the interaction between the atoms to $a_{\mathrm{dc}}$ close to zero
via the wide $s$-wave Feshbach resonance [23], with zero-crossing at 17.1 G.
This way we prepare a BEC soliton. For 4000 atoms the BEC forms a soliton at
interaction $a_{\mathrm{dc}}=-2.7a_{0}$ [24]. To induce the emission of
matter-wave jets, we modulate the magnetic field and therefore the interaction
between the atoms as $a(t)=a_{\mathrm{dc}}+a_{\mathrm{ac}}\sin(\omega t)$ for
time $t_{\mathrm{p}}$, where $a_{\mathrm{ac}}$ and $\omega$ are the amplitude
and the frequency of the modulation, respectively [25]. We let the system
evolve for $t_{\mathrm{e}}$ in the channel and then we take an absorption
picture after a 15 ms time-of-flight. The results are shown in Fig. 1(a). In
addition to previously observed jets of first (J1) and second order (J2) [23],
we observe two pairs of jets with dimensionless velocities

$\displaystyle
u_{\mathrm{G2}}=\frac{v_{\mathrm{G2}}}{v_{\mathrm{0}}}=\pm\frac{1\pm\sqrt{5}}{2},$
(1)

where $v_{\mathrm{0}}=\sqrt{\hbar\omega/m}$ is the velocity of the first-order
jets (J1), $\hbar$ is the reduced Planck constant, and $m$ is the mass of an
atom. We name these “golden jets” (G2) for their golden ratio dimensionless
velocities. The irrational velocity ratio between first order and golden jets
implies the existence of two incommensurable density waves inside the
condensate prior to jet emission. Fig. 1(a) shows three experimental
realizations of golden jets. As can be seen in Fig. 1(a) the golden jets are
not necessarily symmetric around J0 and change in each experimental run
despite the symmetry of involved processes. We attribute this to the asymmetry
in the initial BEC which fluctuates between experimental runs

Table 1:  Jet orders, abbreviations, and formation processes. Beyond $2^{\mathrm{nd}}$ order, only jets observed in the simulations are included.  Order |  Abbrev. |  Formation process
---|---|---
$(u_{1},u_{2})\rightarrow(u^{\prime}_{1},u^{\prime}_{2})$
$0^{\mathrm{th}}$ | J0 | /
$1^{\mathrm{st}}$ | J1 | $(0,0)\rightarrow(-1,1)$
$2^{\mathrm{nd}}$ | J2 | $\pm(1,1)\rightarrow(0,\pm 2)$
G2 | $(0,\pm 1)\rightarrow\pm\left(\frac{1-\sqrt{5}}{2},\frac{1+\sqrt{5}}{2}\right)$
D2 | $(-1,1)\rightarrow\left(-\sqrt{2},\sqrt{2}\right)$
$3^{\mathrm{rd}}$ | J3 | $\pm(2,2)\rightarrow\pm(1,3)$
G3 | $\pm(1,2)\rightarrow\pm\left(\frac{3-\sqrt{5}}{2},\frac{3+\sqrt{5}}{2}\right)$

The origin of golden jets can be understood from the conservation of energy
$\frac{mv_{1}^{2}}{2}+\frac{mv_{2}^{2}}{2}+\hbar\omega=\frac{mv_{1}^{\prime
2}}{2}+\frac{mv_{2}^{\prime 2}}{2}$, and momentum:
$mv_{1}+mv_{2}=mv_{1}^{\prime}+mv_{2}^{\prime}$, where $v_{i}$
($v_{i}^{\prime}$) are the initial (final) velocities of the atoms. A pair of
atoms in J1 is formed from a pair of colliding atoms from the original
condensate (J0) by absorbing one quantum of energy from the modulating
magnetic field. Similarly, a pair J0-J2 forms from a pair of J1 atoms with the
same velocity and a pair G2-G2 from a pair J0-J1. In general the final
velocities of the resulting jets are

$\displaystyle
u^{\prime}_{1,2}=\frac{(u_{1}+u_{2})\pm\sqrt{(u_{1}-u_{2})^{2}+4}}{2},$ (2)

where $u_{1,2}$ are the initial velocities of the atom pair. With each
absorbed photon the number of possible jets increases. By the number of
photons that an atom in a jet has absorbed we define the order of the jets:
the initial BEC is zeroth order, J1 are first order, J2, G2 and others that
form from them are second order and so on, see Fig. 1(b) and Table 1. The
total number of jets increases as a double exponential with order, going as 1,
3, 11, 109, 11605,$\ldots$, and can be approximately calculated as
$1.7949^{2^{n}}$ (from the fit up to $n=4$, where $n$ is the order).
Experimentally, we observe J1, J2 and G2, but we haven’t been able to observe
jets with dimensionless velocity $\pm\sqrt{2}$ (D2, the other second order
jets) or any third order jets. J2 are formed from J1 which move with the same
velocity so the process is equivalent to the primary process in a moving
frame, while the atoms that form G2 move relative to each other, which reduces
the interaction time, making golden jets rare. For D2 the velocity difference
is two times higher, explaining why we don’t observe them.

Figure 2:  (a) Final profile ($t=40$) of the simulation for $\omega=24$,
$k_{\mathrm{dc}}=-0.045$, $k_{\mathrm{ac}}=2.34$, $t_{\mathrm{p}}=8\pi$ and
the fraction of atoms in different jets (calculated as integrals on the
interval (-20,20) around the selected peaks in (b) linear and (c) logarithmic
scale. Bare lines without symbols indicate the supercontinuum regime.  Figure
3:  (a) Phase diagram of jets from GPE simulations. The narrow band between
“no jets” and “J1 + J2” is the region with only J1 jets. (b) Final density
profiles in momentum space $|\psi(q)|^{2}$ ($t=30$) for $\omega=12$ [along the
dotted line in (a)]. (c) and (e) Experimental phase diagram of jets, measured
at $t_{p}=$ 30 ms and 10 ms. The grey lines indicate where the number of atoms
falls by 25% and 50% from the initial number of atoms in the condensate. The
dashed black line in (c) is a fit to the boundary between the ”no jets” and
”J1” region, proportional to $\omega^{1/2}$. The errorbars correspond to the
variance of the atom number and how precisely we can determine the amplitude
corresponding to the boundary between regions. (d) and (f) Averages of 50
absorption images, as examples of different points in the phase diagram
(marked with black stars in the diagrams above). From bottom to top the images
represent regions with: no jets, only J1, J1 and J2, and J1, J2 and G2 (top
two images). The images were taken for evolution time $t_{e}=$ 42.7 ms.
$k_{\mathrm{ac}}=Na_{\mathrm{ac}}/a_{r}$, where $N=4000\pm 200$ is the number
of atoms in the BEC, before we start modulation. $a_{\mathrm{ac}}$ is
calculated from the amplitude of the pulse by a calibration determined using
the Feshbach resonance (details in Supplemental Material).

We turn to numerical simulations of 1D Gross-Pitaevskii equation (GPE) to
better understand the conditions for the formation of each order of jets. Like
in our previous work [23], we start with a soliton as the initial state and
modulate the interaction parameter
$k(t)=Na(t)/a_{r}=k_{\mathrm{dc}}+k_{\mathrm{ac}}\sin(\omega t)$ for a time
$t_{\mathrm{p}}=8\pi$ with different frequencies $\omega$ and amplitudes
$k_{\mathrm{ac}}$. We leave the wavefunction to evolve for $t_{\mathrm{e}}$
and look at the density profile $n(x)=|\psi(x)|^{2}$. In the simulations, time
$t$ is given in units $1/\omega_{r}$, $\omega$ in units $\omega_{r}$, and
distances in units $a_{r}=\sqrt{\hbar/m\omega_{r}}$.

An example of the final density profile can be seen in Fig. 2(a) where all
visible jets are marked. In contrast with the experiment, we observe
additional jets of third order J3 and G3 which form from J2-J2 and J1-J2
pairs, respectively. Surprisingly, G2 and G3 have comparable amplitude even
though G2 jets are second order whereas G3 only arise in third order; this
seems to be due to the large amplitude of J2 jets (from which G3 are formed)
in the simulation.

To quantify the size of jets, we calculate the integrals over a fixed interval
around each peak in the density profile, which gives the fraction of the atoms
in jets of each order (total integral is normalized, $\int|\psi(x)|^{2}dx=1$).
The fraction in the central peak J0 only starts decreasing after the threshold
amplitude for J1 formation is reached, as can be seen in Fig. 2(b). The J1
first increases steeply with amplitude and then levels out due to J0 depletion
and losses to the higher-order processes. It can be seen that only the first
order has a sharp threshold (Fig. 2(b)), while the higher orders increase
smoothly from zero once the first order exists (Fig. 2(c)). D2 are not
observed even in simulations due to the onset of a wide background at high
amplitudes needed to create them (discussed below).

Based on calculated atom fractions, we construct a phase diagram of visibility
of different set of jets, seen in Fig. 3(a). We declare that a new region
begins when the fraction for the next jet exceeds 0.002. We choose this value
so that it is smaller than the fraction of the jets at high modulation
amplitudes and much larger than the noise in the simulation. As the amplitude
$k_{\mathrm{ac}}$ increases, additional jets are enabled and we transition
into new regions. As discussed above the transitions are actually continuous
and not discrete, except from the “no jets” to the “J1” region.

Fig. 3(b) shows the final profiles for a fixed frequency $\omega=12$. The
onset of J1, J2 and other jets can be clearly seen as the modulation amplitude
increases. For very high amplitudes the discrete jets start disappearing in
the increasingly uniform background of a multitude of higher order jets. The
observed spread from the single momentum mode to a wide momentum distribution
resembles the generation of supercontinua in nonlinear optical media, where
various instabilities cause an initially narrow frequency spectrum to disperse
into a wide continuum [26]. The boundary for the supercontinuum can be
determined from the median of the final density profile $n(x)$. It stays close
to zero when we only have discrete jets and when it increases beyond 0.00025
we say that it has transitioned into the supercontinuum region (of course the
transition is actually gradual, the threshold value was chosen to be above
numerical noise and visually without separate jets). The threshold boundary
for J1 formation is proportional to $\omega^{1/2}$, as previously reported
[23, 25]. The fits to lower boundaries for J2 and J3 regions also give
$\omega^{1/2}$ dependences, while the boundaries for golden jets (G2+G3) and
the supercontinuum region go as $\sim\omega^{0.8}$.

We explore the phase diagram of jet formation experimentally as well. Fig.
3(c) and (e) show the experimental phase diagrams of jet formation for
modulation times $t_{p}=$ 10 ms and 30 ms, respectively. We determined the
regions of the phase diagram where we observe J1, J1 and J2, and J1, J2 and
G2.

In contrast to simulation, in the experiments, the number of atoms in jets
varies from shot to shot; for example, see Fig. 1(a). This was previously
observed in 2D [27], where the population in angular modes is also
distributed. More information on the statistics of the number of atoms in each
type of jet is available in the Supplemental Material (Fig. S3). Because of
this variance, gathering statistics is necessary to determine the presence of
the phase at each point in the phase diagram. We choose the boundary to be at
the point where at least one image out of ten shows the corresponding jets.
Examples of averages of 50 absorption images of jets from different parts of
the phase diagram are shown below in Fig. 3(d) and (f) and in the Supplemental
Material (Fig. S2). Note that the distance between the jets and the initial
condensate varies for different modulation amplitudes. This is because jets
are emitted sooner, if the modulation amplitude is higher.

The range of frequencies in the experimental and the theoretical phase diagram
are similar (the frequencies shown in Fig. 3(a) range from 0 to 3.6 kHz,
whereas the experimental phase diagrams in Fig. 3(c) and (e) were measured for
frequencies between 1 kHz and 7 kHz). On the other hand, the values of
$k_{\mathrm{ac}}$ required for excitation of jets are approximately ten times
higher in the experiment than in simulation. This is most probably because the
simulation describes a 1D system, and in the experiment, the BEC is a 3D
object with jets confined into quasi-1D system by a dipole trap beam.

For $t_{p}=$ 30 ms, we are not even able to observe golden jets for modulation
frequencies higher than 3 kHz even though, they are observable using the
shorter $t_{p}=$ 10 ms. In the simulation, the modulation time is
$t_{p}=8\pi$, which corresponds to 44 ms. In the experiment, we used shorter
modulation times, because longer modulation times cause more atom losses and
increase the incoherent background, formed by the excited atoms that are not
part of the jets [23], making the occurrence of jets rarer and more difficult
to observe (Supplemental Material Fig. S1 and S2). However, we see from the
experimental phase diagrams that the shape of the threshold does depend on
$t_{p}$. For a short $t_{p}$ there is no strong dependence on the modulation
frequency, whereas for the longer $t_{p}$ the threshold boundary for the
formation of jets is proportional to $\omega^{1/2}$, as expected from the
simulation and previous experiments [23]. Additionally, the longer modulation
times sometimes cause the jets to be emitted twice, making the interpretation
of the absorption images less reliable.

As we increase the amplitude in the experiment, the total number of atoms
starts decreasing and the incoherent background between the jets becomes
stronger, as seen in Fig. 3, where we also plot the thresholds at which the
total number of atoms is 25% and 50% smaller than the number of atoms in the
initial BEC. The decrease in the number of atoms is likely due to partial
collapses of the BEC, caused by the strong attractive interactions the BEC is
subjected to during the modulation. The momentum distribution of the
incoherent background also remains of the same width for increasing amplitude
and does not increase beyond the position of J2 jets as the supercontinuum
would. Atom losses and increased fraction of incoherent background, both
absent from the simulation, decrease the effective modulation amplitude
preventing the supercontinuum from forming.

In this work, we show the formation of matter-wave jets with golden ratio
dimensionless velocity $u_{\mathrm{G2}}={(1+\sqrt{5})/2}$, which form when a
BEC in a quasi-1D geometry is subjected to a single-frequency interaction
modulation. These jet velocities and the correspoding density waves are not
commensurate with previously observed jets with dimensionless velocities
$u_{\mathrm{J1}}=1$ and $u_{\mathrm{J2}}=2$. It needs to be noted, that the
occurrence of incommensurate density waves can easily be achieved with two-
frequency driving, but here it emerges from single-frequency driving.
Numerically and experimentally, we explore the formation of incommensurable
jets for a wide range of modulation frequencies and amplitudes, resulting in a
phase diagram of possible jets, while also uncovering an unexpected
supercontinuum region. The supercontinuum likely arises due to the double-
exponential growth of the number of possible modes, while the wavenumber only
grows linearly, resulting in a densely populated momentum space. The
supercontinuum is potentially experimentally observable if the three-body
losses and the incoherent background could be reduced by increasing the radial
trapping frequency, making the system more 1D. The demonstrated formation of
incommensurable matter-wave jets could, for example, be used in multispecies
experiments. The incommensurable density waves inside one species could be
used as a disordered potential for the other species, if the jets are
prevented from escaping the trap (provided there exist appropriate Feshbach
resonances). This would eliminate the need for an external disordered
potential and even a regular lattice that is usually needed to capture
impurities, which act as a source of disorder [12, 13, 14, 15, 16, 17, 18, 19,
20].

###### Acknowledgements.

We thank Cheng Chin and Matjaž Gomilšek for helpful discussions. We would also
like to thank Samo Beguš and Davorin Kotnik for their help with electronics.
This work was supported by the Slovenian Research Agency (research core Grants
No. P1-0125, No. P1-0099 and P1-0416, and research projects No. J2-8191 and
No. J2-2514).

## References

  * Bak [1982] P. Bak, Reports on Progress in Physics 45, 587 (1982).
  * Lifshitz [2007] R. Lifshitz, Zeitschrift für Kristallographie 222, 313 (2007).
  * Janssen [2012] T. Janssen, Acta Crystallographica Section A Foundations of Crystallography 68, 667 (2012).
  * Feng  _et al._ [2013] Y. Feng, J. Wang, D. M. Silevitch, B. Mihaila, J. W. Kim, J.-Q. Yan, R. K. Schulze, N. Woo, A. Palmer, Y. Ren, J. van Wezel, P. B. Littlewood, and T. F. Rosenbaum, Proceedings of the National Academy of Sciences 110, 3287 (2013).
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Wang  _et al._ [2018] Y. Wang, B. Ge, Y. Cui, R. Yu, S. Chen, and Z. Ren, Materials Today Physics 5, 7 (2018).
  * Miao  _et al._ [2019] H. Miao, R. Fumagalli, M. Rossi, J. Lorenzana, G. Seibold, F. Yakhou-Harris, K. Kummer, N. B. Brookes, G. D. Gu, L. Braicovich, G. Ghiringhelli, and M. P. M. Dean, Physical Review X 9, 031042 (2019).
  * Storeck  _et al._ [2020] G. Storeck, J. G. Horstmann, T. Diekmann, S. Vogelgesang, G. von Witte, S. V. Yalunin, K. Rossnagel, and C. Ropers, Structural Dynamics 7, 034304 (2020).
  * Carr  _et al._ [2020] S. Carr, S. Fang, and E. Kaxiras, Nature Reviews Materials 5, 748 (2020).
  * Lu  _et al._ [2021] X. Lu, B. Lian, G. Chaudhary, B. A. Piot, G. Romagnoli, K. Watanabe, T. Taniguchi, M. Poggio, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, Proceedings of the National Academy of Sciences 118, 10.1073/pnas.2100006118 (2021).
  * Hu  _et al._ [2022] Y. Hu, T. Zhang, D. Zhao, C. Chen, S. Ding, W. Yang, X. Wang, C. Li, H. Wang, D. Feng, and T. Zhang, Nature Communications 13, 10.1038/s41467-022-28104-2 (2022).
  * Lye  _et al._ [2007] J. E. Lye, L. Fallani, C. Fort, V. Guarrera, M. Modugno, D. S. Wiersma, and M. Inguscio, Physical Review A 75, 061603(R) (2007).
  * Roati  _et al._ [2008] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, Nature 453, 895 (2008).
  * Fallani  _et al._ [2008] L. Fallani, C. Fort, and M. Inguscio, in  _Advances In Atomic, Molecular, and Optical Physics_ (Elsevier, 2008) pp. 119–160.
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider, and I. Bloch, Physical Review X 7, 041047 (2017).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Science 349, 842 (2015).
  * D’Errico  _et al._ [2014] C. D’Errico, E. Lucioni, L. Tanzi, L. Gori, G. Roux, I. P. McCulloch, T. Giamarchi, M. Inguscio, and G. Modugno, Physical Review Letters 113, 095301 (2014).
  * Gadway  _et al._ [2011] B. Gadway, D. Pertot, J. Reeves, M. Vogt, and D. Schneble, Physical Review Letters 107, 145306 (2011).
  * Lin  _et al._ [2019] J. Lin, J. Nan, Y. Luo, X.-C. Yao, and X. Li, Physical Review Letters 123, 233603 (2019).
  * Reeves  _et al._ [2014] J. B. Reeves, B. Gadway, T. Bergeman, I. Danshita, and D. Schneble, New Journal of Physics 16, 065011 (2014).
  * Fu  _et al._ [2018] H. Fu, L. Feng, B. M. Anderson, L. W. Clark, J. Hu, J. W. Andrade, C. Chin, and K. Levin, Phys. Rev. Lett. 121, 243001 (2018).
  * Zhang  _et al._ [2020] Z. Zhang, K.-X. Yao, L. Feng, J. Hu, and C. Chin, Nat. Phys. 16, 652 (2020).
  * Mežnaršič  _et al._ [2020] T. Mežnaršič, R. Žitko, T. Arh, K. Gosar, E. Zupanič, and P. Jeglič, Phys. Rev. A 101, 031601(R) (2020).
  * Mežnaršič  _et al._ [2019] T. Mežnaršič, T. Arh, J. Brence, J. Pišljar, K. Gosar, Ž. Gosar, R. Žitko, E. Zupanič, and P. Jeglič, Phys. Rev. A 99, 033625 (2019).
  * Clark  _et al._ [2017] L. W. Clark, A. Gaj, L. Feng, and C. Chin, Nature 551, 356 (2017).
  * Dudley  _et al._ [2006] J. M. Dudley, G. Genty, and S. Coen, Reviews of Modern Physics 78, 1135 (2006).
  * Hu  _et al._ [2019] J. Hu, L. Feng, Z. Zhang, and C. Chin, Nat. Phys. 15, 785 (2019).