# Self-supervised Multidimensional Scaling with $F$-ratio: Improving
Microbiome Visualization

Hyungseok Kim∗ Massachusetts Institute of Technology  Soobin Kim∗ University
of California, Davis  Megan M. Morris  Lawrence Livermore National Laboratory
Jeffrey A. Kimbrel  Lawrence Livermore National Laboratory  Xavier Mayali
Lawrence Livermore National Laboratory  Cullen R. Buie† Massachusetts
Institute of Technology

###### Abstract

Multidimensional scaling (MDS) is an unsupervised learning technique that
preserves pairwise distances between observations and is commonly used for
analyzing multivariate biological datasets. Recent advances in MDS have
achieved successful classification results, but the configurations heavily
depend on the choice of hyperparameters, limiting its broader application.
Here, we present a self-supervised MDS approach informed by the dispersions of
observations that share a common binary label ($F$-ratio). Our visualization
accurately configures the $F$-ratio while consistently preserving the global
structure with a low data distortion compared to existing dimensionality
reduction tools. Using an algal microbiome dataset, we show that this new
method better illustrates the community’s response to the host, suggesting its
potential impact on microbiology and ecology data analysis 1.

Keywords: Multidimensional scaling, microbiome, F-statistic, dimensionality
reduction.

**footnotetext: Equal contribution.$\dagger$$\dagger$footnotetext:
Correspondence<EMAIL_ADDRESS>Code available at
https://github.com/hyu-kim/mds-hypothesis-testing.

##  1 Introduction

Biotechnological advances in the past several decades have expanded the size
and features of the multivariate data, necessitating a dimensionality
reduction as a tool for the interpretation. By extracting an essential
information from the biological data, the dimensionality reduction seeks a
visual representation of the multivariate in a lower dimensional space. To
retain a consistent data structure while performing the dimension reduction, a
configuration is sought in a way that preserves dissimilarity between samples,
a process called as the multidimensional scaling (MDS). Compared to other
nonlinear methods in dimensionality reduction [21, 18, 22], MDS is known to
retain a global structure and represent a long-range interaction between
samples, allowing its popularity for a long time since its inception [8].

In the classical MDS, the configuration is determined in a way that minimizes
the difference between the dissimilarity in the original and the lower
dimensional space. Because biological observations often present a
compositional structure (highly skewed, zero-inflated [11]), choosing an
appropriate dissimilarity metric is an important step for processing the data
and representing in the low dimensional space. For example, when interpreting
a microbial community dataset such as 16S rRNA gene expression, a distance
metric such as Unifrac [16, 17] allows to incorporate taxonomic information as
well as the compositional structure.

Whereas the dimensionality reduction provides a qualitative interpretation on
the multivariate data structure, statistical inference such as hypothesis
testing is carried out at the same time to give a quantitative analysis on the
structural difference by biological responses. One way to test the difference
between response groups is to compare dispersions of the multivariate
structure that are calculated across or within the sample groups, which gives
an F-statistic[11]. The process employs a different information than those
used for MDS, in that it is unsupervised learning method and does not use
label information. Therefore, MDS configuration under a low dimension (e.g.,
2) may not fully explain a small but statistically meaningful difference
between groups of different treatment.

The insufficient explanation by the classical MDS encourages to revise the
method and to address a structural hypothesis by including an external
information which is conferred by the responses or class labels [10]. Broadly
termed as the confirmatory MDS, this approach is based on applying an external
constraint to the classical MDS, thereby providing a more contextual
illustration of data structure [2]. Given that several local configurations
can be produced from classical MDS, there is a notion accepting altered
configurations up to a point where they do not deviate too much from the
original configuration.

In these confirmatory MDS methods, an objective function is constructed by
adding a confirmatory term to the stress, weighed by a hyperparameter, and
informed by the labels or responses. These confirmatory MDS methods have
successfully visualized the multivariate structure in a way that
differentiates sample groups with a discriminative purpose [6, 24]. However,
choosing a proper hyperparameter remains as a bottleneck towards a broader
application of these inventive methods. For example, setting a high
hyperparameter can result in an undesirable stress and a misleading
configuration heavily distorted from the original.

In this work we introduce a weakly supervised version of MDS informed by
F-statistic of multivariate dataset, which we refer to as FMDS. In other
words, FMDS is informed by a hypothesis testing results under a binary class
setting using $F$-statistic. The proposed approach is motivated by a purpose
to explain a statistical difference between groups, if any, and to incorporate
it to the confirmatory analysis by combining with classical MDS configuration.
Our motivation distinguishes itself from existing variations in confirmatory
MDS because the method does not target to directly discriminate between
groups. Furthermore, by characterizing the proposed framework we show that
FMDS configuration is less dependent to the choice of the model
hyperparameter, mediating the existing issue with the distortion.

##  2 Related Work

##### Classical multidimensional scaling:

Consider a balanced design where the number of total observations is $N$, and
each observation $x_{i}$ is $S$-dimensional, pertaining to a set of labels
$y_{i}\in\\{0,1\\}$ for every $i=1,\cdots N$. Using on the set of observations
$(x_{1},\cdots x_{N})$, a pairwise distance $d_{ij}\in\\{0,\infty\\}$ between
$x_{i}$ and $x_{j}$ can be obtained using an appropriate distance metric
(e.g., Euclidean, Bray-Curtis [4], UniFrac [16]). In the classical MDS, a
lower-dimensional configuration
$\mathbf{z}=(\mathbf{z}_{1},\cdots\mathbf{z}_{N})\in\operatorname*{\mathbb{R}}^{N\times
2}$ is sought in a way that preserves the pairwise distance while the
dimension is reduced. This is enabled by minimizing the “raw stress” [3]:
${\rm
Stress}=\frac{1}{2}\sum_{i,j}(d_{ij}-\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2})^{2}$.
Therefore, the classical MDS is unsupervised learning and it does not require
a set of labels $y_{i}$ for training.

##### Supervised multidimensional scaling [23]:

Supervised multidimensional scaling imposes an additional constraint to the
configuration using a confirmatory feature that is informed by class labels.
The purpose is to discriminate observations by the labels as well as carrying
out the task of minimizing the stress. To achieve this an objective function
$O(\mathbf{z})$ is proposed by summing the raw stress and a confirmatory term,

$O(\mathbf{z})=\frac{1}{2}(1-\alpha)\sum_{i,j}(d_{ij}-\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2})^{2}+\alpha\sum_{i,j:y_{j}>y_{i}}(y_{j}-y_{i})\sum_{s=1}^{2}\left(\frac{D_{ij}}{\sqrt{2}}-(z_{js}-z_{is})\right)^{2},$
(1)

where the confirmatory term contains labels $y_{i}$ so that the configuration
points $\mathbf{z}_{i}$ and $\mathbf{z}_{j}$ are not too far apart when
$y_{i}\approx y_{j}$. In addition, a hyperparameter $\alpha$ is chosen to
determine the degree of classification task over the stress minimization.

##### Non-parametric multivariate analysis of variance:

When testing a difference between groups of multivariate data, test statistics
are obtained through calculating variances or a sum of distances /
dissimilarities. Most widely used metric is $F$-statistic, a ratio between two
sums respectively calculated from inter- and intra-group variances. However,
performing $F$-test necessitates an assumption that the observations follow a
normal distribution with common variance, which is hardly met in practice when
dealing with biological dataset. In this non-parametric context, an
alternative approach has been proposed by permuting the labels iteratively to
create an empirical distribution [1, 12]. First, define the pseudo $F$-ratio
as

$F=\frac{\sum_{i,j}d_{ij}^{2}-2\sum_{i,j}{\rm 1}\kern-2.40005pt{\rm
I}\\{y_{i}=y_{j}\\}d_{ij}^{2}}{2\sum_{i,j}{\rm 1}\kern-2.40005pt{\rm
I}\\{y_{i}=y_{j}\\}d_{ij}^{2}}\cdot(N-2),$ (2)

where ${\rm 1}\kern-2.40005pt{\rm I}\\{\cdot\\}$ denotes an indicator
function. While the pseudo $F$ does not always follow the $F$-distribution
under the non-parametric setting, an empirical distribution can be constructed
instead by ‘permuting’ the labels for a large enough size of dataset [11, 1].
In other words, denote $F^{\pi}$ as a new $F$-ratio that is obtained from the
permuted labels, and by iterating the permutation to obtain a $P$-value,

$P=\frac{\textrm{Number of case where }(F^{\pi}\geq F)}{\textrm{Number of
total repeat}}.$ (3)

The procedure was first proposed by Anderson [1], known as the permutational
multivariate analysis of variance (PERMANOVA).

##  3 Approach

We have set our goal to incorporate the multivariate hypothesis testing result
for determining an MDS configuration of the dataset. To achieve the goal, we
take a self-supervised approach with the following two steps. First, a two-
dimensional configuration is trained by performing the classical MDS (i.e.,
unsupervised). Next, the configuration points are adjusted in a way that best
represents a $P$ value calculated from the multivariate, binary-labeled
dataset. This is enabled by adding a confirmatory term to the raw stress, a
similar approach proposed by the supervised MDS [23],

$O(\mathbf{z})=\underbrace{\frac{1}{2}\sum_{i,j}(d_{ij}-\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2})^{2}}_{\text{raw
stress}}+\lambda\cdot\underbrace{\frac{1}{2}\left|\sum_{i,j}\left[1-2\epsilon_{ij}\left(1+\frac{f_{\mathbf{z}}(F)}{N-2}\right)\right]\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2}^{2}\right|}_{\text{{F}-informed
confirmatory term}},$ (4)

where
$f_{\mathbf{z}}(F):\operatorname*{\mathbb{R}}\rightarrow\operatorname*{\mathbb{R}}$
is a scalar function that maps $F$-ratio from $S$-dimension onto the lower
dimension, 2, which is determined by the configuration $\mathbf{z}$.
Introducing the mapping $f_{\mathbf{z}}$ allows to address a difference in
$F$-ratio distributions between $S$\- and 2-dimension, thereby bridging a gap
in the P values. A detailed derivation of $f_{\mathbf{z}}$ and a description
on the confirmatory term is provided in Supplementary Note 1, and its
pseudocode is shown in Algorithm 1.

1 Function  _mapping $(\mathbf{z},d,y)$_:

Input: $\mathbf{z}$: | two-dimensional configuration
---|---
$d$: | pairwise distance matrix
$y$: | labels set

Output: $f_{\mathbf{z}}$: | mapping function
---|---
$L$: | sorted list of permuted $F$
$L_{\mathbf{z}}$: | sorted list of permuted $F_{\mathbf{z}}$

2 for  _$1\leq i\leq 999$_ do

3 $y^{\pi}_{1}\leftarrow$ shuffle $y$ by a random permutation

4 $F^{\pi}\leftarrow$ $F$-ratio by $d,y^{\pi}_{1}$ (Equation S1)

5 $L\leftarrow$ append $F^{\pi}$ to $L$

6 $y^{\pi}_{2}\leftarrow$ shuffle $y$ by another random permutation

7 $F^{\pi}_{\mathbf{z}}\leftarrow$ $F$-ratio by $\mathbf{z},y^{\pi}_{2}$
(Equation S2)

8 $L_{\mathbf{z}}\leftarrow$ append $F^{\pi}_{\mathbf{z}}$ to $L_{\mathbf{z}}$

9

10 end for

11 $L\leftarrow\text{sort}(L)$

12 $L_{\mathbf{z}}\leftarrow\text{sort}(L_{\mathbf{z}})$

13 $f_{\mathbf{z}}\leftarrow\text{LOESS}(L,L_{\mathbf{z}})$

14 return $f_{\mathbf{z}}$

Algorithm 1 Mapping from $F$ to $F_{\mathbf{z}}$ with random permutation.

###  3.1 Majorize-Minimization (MM) algorithm

Minimizing raw stress to perform MDS has been achieved by several iteration
methods in the past few decades, part of which have originated in the field of
graph drawing [13]. Because the stress is a non-convex function in terms of
$\mathbf{z}$, there is no global minimum and it may require a different
initialization [8, 26]. This has also led with proposing several optimization
methods including 2D Newton–Raphson [13], (stochastic [26]) gradient descent
[15], divide-and-conquer [25, 19], and majorization [7].

In our problem setting, we implemented the majorization (or Majorize-
Minimization) approach to minimize the FMDS objective (Equation 4). In brief,
for every $k=1,\cdots N$, Equation 4 can be minimized by writing

$\displaystyle\mathbf{z}_{k}^{*}$
$\displaystyle=\operatorname*{arg\,min}_{\mathbf{z}_{k}}O(\mathbf{z})$ (5)
$\displaystyle=\operatorname*{arg\,min}_{\mathbf{z}_{k}}\,\sum_{j=1}^{N}\left[1+\lambda\delta(\mathbf{z})\left(1-2\epsilon_{jk}\left(1+\frac{f_{\mathbf{z}}(F)}{N-2}\right)\right)\right]\|\mathbf{z}_{k}-\mathbf{z}_{j}\|_{2}^{2}-2d_{jk}\|\mathbf{z}_{k}-\mathbf{z}_{j}\|_{2}$
(6)

where we define $\epsilon_{ij}$ and $\delta(\mathbf{z})$ as

$\displaystyle\epsilon_{ij}$ $\displaystyle={\rm 1}\kern-2.40005pt{\rm
I}\\{y_{i}=y_{j}\\},\quad\delta(\mathbf{z})=\text{sign}\left\\{\,\sum_{i,j}\left[1-2\epsilon_{ij}\left(1+\frac{f_{\mathbf{z}}(F)}{N-2}\right)\right]\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2}^{2}\right\\}.$
(7)

Majorizing with Equation 6 gives a quadratic expression in terms of
$\mathbf{z}$, as similarly described by Borg [3] and Witten [23]. Finally, a
derivative is analytically taken with respect to $z_{ks}$ and is set to zero
to find a local minimum, given that the expression is convex with assumptions.
Detailed procedure is described in Supplementary note 2 and we provide the
update rule in Algorithm 2.

1 Function  _mmfmds $(\lambda,d,y)$_:

Input: $\lambda$: | hyperparameter balancing raw stress and confirmatory term
---|---
$d$: | pairwise distance matrix
$y$: | labels

Output: $\mathbf{z}$: | two-dimensional configuration
---|---

2 for  _$1\leq t\leq T$_ do

3 for  _$1\leq k\leq N$_ do

4 Apply Equation 2 and compute $F$.

5 Apply Equation 7 and compute $\delta(\mathbf{z}).$

6 Apply the mapping algorithm (1) and compute $f_{\mathbf{z}}(F)$.

7
$\displaystyle\mathbf{z}_{k}\leftarrow\frac{(N-2)}{N(N-2)-N\lambda\delta(\mathbf{z})f_{\mathbf{z}}(F)}\times$
$\displaystyle\qquad\left\\{\sum_{j=1}^{N}\,\left[1+\lambda\delta(\mathbf{z})\left(1-2\epsilon_{jk}\left(1+\frac{f_{\mathbf{z}}(F)}{N-2}\right)\right)\right]\mathbf{z}_{j}+d_{jk}\frac{\mathbf{z}_{k}-\mathbf{z}_{j}}{\|{\mathbf{z}}_{k}-\mathbf{z}_{j}\|_{2}}\right\\}$

8 end for

9

10 end for

11

12 return $\mathbf{z}$

Algorithm 2 Majorize-Minimization for FMDS

##  4 Experiments

We first characterized the behavior of FMDS across a range of hyperparameters
and compared its performance to benchmark dimensionality reduction methods,
including supervised MDS [23], UMAP [18], and a self-supervised neural network
such as SimCLR [5]. For the evaluation we considered two types of datasets,
simulated and experimental, both of which serve as examples where the two-
dimensional configurations from classical MDS were not consistent with the
hypothesis testing results. Next, in Section 4.2, we demonstrate how FMDS can
improve visualizing biological samples in a way that multivariate $F$ test
results are addressed.

###  4.1 Evaluating performance of F-informed MDS

We evaluated the performance of FMDS by quantifying a degree of deviation of
the configuration $\mathbf{z}$ from the original distances $d_{ij}$. It is
carried out by calculating two metrics; the first is the “normalized” stress
or Stress-1 [3],

$\text{Stress-1}=\frac{\sum_{i,j}(d_{ij}-\|\mathbf{z}_{i}-\mathbf{z}_{j}\|_{2})^{2}}{\sum_{i,j}d_{ij}^{2}},$
(8)

and the second using a Shepard diagram [9] and its correlation coefficient.

####  4.1.1 Simulated dataset

We first considered a three-dimensional dataset of two balanced groups that
follow normal distributions slightly different in means, $\mu_{0}$, $\mu_{1}$,
but the same covariance matrix $D$. Such set of observations $x_{i}$ can be
constructed as, for example,

$x_{i}\sim\begin{cases}\mathcal{N}\left(\mu_{0},D\right),&i=1,2,\cdots 50\\\
\mathcal{N}\left(\mu_{1},D\right),&i=51,52,\cdots 100,\end{cases}$ (9)

where $\mu_{0}=[0,0,0]^{\top}$, $\mu_{1}=[0,0,1]^{\top}$, and
$D=[[3,0,0],[0,3,0],[0,0,1]]$.

Here the means are different at the third dimension, where the lowest variance
was imposed among the principal diagonals of $D$ (Supplementary Figure 1).
Therefore, a classical MDS does not distinguish groups in a two-dimensional
configuration, but performing a multivariate hypothesis testing on the
datasets, e.g., PERMANOVA [1], indicates there is a statistically significant
difference between the groups with $P=0.005$.

Using the simulated dataset, we measured the normalized stress (Equation 8)
and correlation in Shepard diagram. Using Algorithm 2 we confirmed the
configuration converged over a number of iterations for a range of
hyperparameter up to $\lambda=0.7$, a value that ensured the confirmatory term
in Equation 4 becomes negligible compared to MDS term (Figure 1).

Figure 1: Comparison of MDS and confirmatory terms by iteration of Algorithm 2
for a range of $\lambda$. Each term values were calculated using simulated
dataset following Equation 9.

Figure 2 shows a summarized result of the performance of FMDS compared to
existing MDS methods. Regardless of a choice of the hyperparameter $\lambda$,
FMDS visualization exhibited a consistent performance with its stress $\sim
0.2$ and correlation coefficient $>0.9$ (Figure 2a). The behavior is in
contrast to the previous supervised MDS [23] which was monotonically dependent
to $\lambda$, as expected, because SMDS distinguishes groups at the expense of
the original distance structure. Overall, the stress obtained from the
proposed MDS is consistently lower than those from SMDS. Similarly, as
displayed in Figure 2b, Shepard plot shows the proposed FMDS presents a higher
correlation of the sample pair distance in between the original and two-
dimensional space.

Figure 2: Evaluation of proposed FMDS using simulated data. (a) Performance of
FMDS compared to supervised MDS [23] by measuring stress and Pearson
correlation coefficient from Shepard plot. (b) Shepard plot of the proposed
MDS comparing to SuperMDS for a hyperparameter $\lambda=0.5$. More Shepard
plots are provided in Supplementary Figure 2.

We then compare the performance of FMDS to other dimensionality reductions
such as UMAP [18], that are widely used in visualizing multivariate data.
While UMAP does not employ the same hyperparameter as FMDS or SMDS, we
evaluated its performance by varying the size of local neighborhood, used for
manifold approximation, from 5 to 30. For a better comparison supervised
learning mode with UMAP is also implemented. Table 1 shows the summarized
result of performance by FMDS, SMDS, and UMAP with two different modes for the
range of hyperparameters using simulated dataset. The result suggests the
proposed FMDS produces a 2D configuration by imposing a consistently lower
distortion of distance structure compared to other methods.

Table 1: Summary of performance of dimensionality reduction methods using simulated data. | FMDS | SMDS [23] | UMAP [18] | UMAP (supervised) [18]
---|---|---|---|---
Stress-1 | 0.15 – 0.20 | 0.13 – 0.42 | 0.32 – 0.44 | 0.61 – 0.70
Pearson correlation | 0.90 – 0.96 | 0.70 – 0.96 | 0.69 – 0.77 | 0.10 – 0.30

####  4.1.2 Microbial community

We next take a microbial community as another dataset to compare the
performances of different dimensionality reductions. The dataset presents a
compositional structure expressed by abundance of 16S rRNA gene of 72
bacterial taxa. As a distance metric, the weighted Unifrac [17] is chosen to
obtain pairwise distance between individual community samples. We considered
two datasets (Site 1 and Site 2) where each contains thirty-six, balanced
community samples with a binary label (e.g., with or without a presence of
bacterial host) [14].

The performance of FMDS is evaluated by calculating Stress-1 and the Shepard
plot visualization. Again we observe that FMDS configuration produced a low
stress that is less dependent on the hyperparameter $\lambda$, which even
decreases with nonzero $\lambda$ (0.1, 0.3) than a classical MDS (Figure 3a).
Shepard plot and Pearson correlation also show that the configurations nicely
preserve the original distance in the microbial community data, except for a
case when the largest $\lambda$ = 0.5 is applied to Site 1 community dataset.

Figure 3: Evaluation of the FMDS using microbial community data, measured by
(a) Stress-1 and Pearson correlation from Shepard plot. (c) Shepard plot from
each sample site. More results are provided in Supplementary Figure 4.

##### Self-supervised learning with SimCLR and PopPhy-CNN.

In addition to the above, we also sought to compare our FMDS with existing
neural network models that are used for dimensionality reduction. For
converting compositional microbial abundance with its phylogenetic information
into a matrix, we implement PopPhy-CNN [20] architecture. In brief, the
encoder consists of one Gaussian noise filter, two 2D convolution layers (with
kernel size of 5 by 3), and one fully connected layer with 32 output nodes.
For the self-supervised learning framework we choose SimCLR [5], where the
data augmentation is performed by applying random brightness and contrast
filter. In a pretraining step of SimCLR, a model is constructed by compiling
encoder, projection head (two dense layers each of 32 output nodes), and one
dense layer (10 output nodes). Site 1 and 2 datasets are merged and split into
train and evaluation data (60, 12 each) for the pretraining step, resulting a
linear probing accuracy of 53.3% after 50 training epochs. In the following
finetuning step, the encoder is added with linear probe, which resulted in a
validation accuracy of 83.3% after 50 epochs. Trained encoder is used to
obtain a 32 nodes-sized feature for each microbial community sample in Site 1
and 2. Detailed parameters and steps for training / evaluating the neural
network architecture is described in Supplementary note 3.

##### Performace evaluation of FMDS with other methods.

Table 2 shows the summarized performance of FMDS with microbiome dataset
compared to existing dimensionality reduction methods (UMAP, neural network
architecture). Again we observe that the proposed MDS delivers a two-
dimensional configuration with less distorted distance structure when compared
to self-supervised learning model (SimCLR) or supervised UMAP. Overall, the
evaluation suggests that FMDS outperforms existing (self-) supervised learning
tools by preserving the global structure in regards to the pairwise distances.

Table 2: Summary of dimensionality reduction tool performances using 72 microbiome samples. |  | FMDS
---
(proposed)
| SimCLR [5] +
---
PopPhy-CNN [20]
UMAP | UMAP (supervised)
Stress-1 | 0.21 – 0.40 | 0.54 | 0.31 – 0.54 | 0.63 – 0.69
Pearson correlation | 0.61 – 0.92 | 0.076 | 0.63 – 0.69 | 0.06 – 0.29

###  4.2 Discriminant analysis with FMDS

We next demonstrate how FMDS can handle two-dimensional configuration by
addressing sample group difference. First, we consider the simulated dataset
(Section 4.1.1) again where a classical MDS does not distinguish the binary
groups in two-dimensional configuration ($P=0.914$) because the group
difference lies in the third dimension which is not identifiable (Figure 4a).
However, the group difference becomes more visible when FMDS is employed, as
shown in Figure 4b,c. The distinctions are also verified by a low $P$-value
resulting from PERMANOVA test using the two-dimensional configurations
($P=0.003$).

Figure 4: 2D visualization of simulated data using MDS. Configurations are
obtained by setting a hyperparameter (a) $\lambda=0$ (classical MDS), (b)
$\lambda=0.3$, and (c) $\lambda=0.5$. For each configuration, a $p$-value is
given based on PERMANOVA test.

##### Bacterial community dataset.

Finally, we demonstrate how FMDS visualization can improve interpreting
microbiome structure while addressing $F$-test results at the same time. Our
datasets, Site 1 and 2, present a unique case where a classical two-
dimensional MDS does not fully explain statistical test results on group
differences. As shown in Figure 5a, groups in site 1 are dispersed in a
different location whereas site 2 groups are not, when visualized using the
classical MDS. In both Sites, however, moderately small $P$-values are
obtained ($<0.1$, Table 3), indicating the group difference in the community
structure is, in fact, statistically significant.

We then visualize the microbiome data using classical and FMDS and compare the
configurations. As expected, for Site 1 community samples the configuration
retains its distinction between the class labels regardless of the choice of
the hyperparameter $\lambda$ (Figure 5). Moreover, for Site 2 samples we
observe a higher distinction between the groups with increasing $\lambda$ when
compared to classical MDS (Figure 5b, Supplementary Figure 3). The observation
with the visualizations is justified by a quantitative measure using $P$-value
calculated on the 2D configurations (Table 3).

Figure 5: Visualization of algal microbiome data using FMDS. (a) Classical and (b) proposed FMDS comparing two-dimensional configurations sampled at Site 1 and 2. More configurations with other value of hyperparameters are displayed in Supplementary Figure 3. Table 3: Statistical significance on the group difference between two treatments using PERMANOVA test [1]. $P$-value | Site 1 | Site 2
---|---|---
Original |  $<$0.001 | 0.093
Classical MDS |  $<$0.001 | 0.05
FMDS ($\lambda=0.5$) |  $<$0.001 | 0.099

##  5 Conclusion

In this work, a self-supervised multidimensional scaling is proposed based on
the F-statistic and is characterized by comparing the configuration to
hypothesis testing results. Using simulated dataset and ecological samples, it
is shown that the proposed MDS outperforms an existing method for addressing
class labels, as evaluated by its stress and Shepard plot. One limitation of
FMDS is that it requires a high computational cost compared to existing
methods, in its algorithmic performance during iterations (Supplementary note
4).

A special attention has been given on its behavior where its performances were
less dependent on the choice of hyperparameter to an extent where algorithm
effectively converges. The finding suggests that the new approach mediates the
downside of typical confirmatory MDS as a dimensionality reduction tool, in
which stress minimization has been underscored by a discriminatory purpose via
a high value of hyperparameter setting. In practice, users may avoid the
hassle of hyperparameter selection using such as cross-validation.

Taken together, The proposed MDS can be useful for interpreting a biological
dataset with more visual information to explain a hypothesis testing result at
the same time, thereby presenting a broader applicability of MDS in modern
biological and multivariate data analysis.

## Acknowledgements

The work was supported by the Department of Energy’s Genome Sciences Program
grant SCW1039. H.K. was partly supported by the Kwanjeong Educational
Foundation. S.K. was partly supported by Global Korea Scholarship (GKS).

## References

  * [1] Anderson, M. J. A new method for non-parametric multivariate analysis of variance.  Austral Ecology 26, 1 (2001), 32–46.
  * [2] Borg, I., and Groenen, P. Confirmatory MDS.  Springer New York, New York, NY, 1997, pp. 181–197.
  * [3] Borg, I., and Groenen, P. MDS Models and Measures of Fit.  Springer New York, New York, NY, 2005, pp. 37–61.
  * [4] Bray, J. R., and Curtis, J. T. An ordination of the upland forest communities of southern wisconsin.  Ecological Monographs 27, 4 (1957), 325–349.
  * [5] Chen, T., Kornblith, S., Norouzi, M., and Hinton, G. A simple framework for contrastive learning of visual representations.  In Proceedings of the 37th International Conference on Machine Learning, vol. 119, Proceedings of Machine Learning Research, pp. 1597–1607.
  * [6] Cox, T. F., and Ferry, G. Discriminant analysis using non-metric multidimensional scaling.  Pattern Recognition 26, 1 (1993), 145–153.
  * [7] de Leeuw, J. Convergence of the majorization method for multidimensional scaling.  Journal of Classification 5, 2 (1988), 163–180.
  * [8] Demaine, E., Hesterberg, A., Koehler, F., Lynch, J., and Urschel, J. Multidimensional scaling: Approximation and complexity, September 01, 2021 2021.
  * [9] Dexter, E., Rollwagen-Bollens, G., and Bollens, S. M. The trouble with stress: A flexible method for the evaluation of nonmetric multidimensional scaling.  Limnology and Oceanography: Methods 16, 7 (2018), 434–443.
  * [10] Ding, C. S. Testing Pattern Hypotheses with MDS.  Springer International Publishing, Cham, 2018, pp. 165–173.
  * [11] Gijbels, I., and Omelka, M. Testing for homogeneity of multivariate dispersions using dissimilarity measures.  Biometrics 69, 1 (2013), 137–145.
  * [12] Holmes, A. P., Blair, R. C., Watson, J. D., and Ford, I. Nonparametric analysis of statistic images from functional mapping experiments.  Journal of Cerebral Blood Flow & Metabolism 16, 1 (1996), 7–22.
  * [13] Kamada, T., and Kawai, S. An algorithm for drawing general undirected graphs.  Information Processing Letters 31, 1 (1989), 7–15.
  * [14] Kim, H., Kimbrel, J. A., Vaiana, C. A., Wollard, J. R., Mayali, X., and Buie, C. R. Bacterial response to spatial gradients of algal-derived nutrients in a porous microplate.  The ISME Journal 16, 4 (2022), 1036–1045.
  * [15] Kruskal, J. B. Nonmetric multidimensional scaling: A numerical method.  Psychometrika 29, 2 (1964), 115–129.
  * [16] Lozupone, C., and Knight, R. Unifrac: a new phylogenetic method for comparing microbial communities.  Applied and Environmental Microbiology 71, 12 (2005), 8228–35.
  * [17] Lozupone, C. A., Hamady, M., Kelley, S. T., and Knight, R. Quantitative and qualitative beta diversity measures lead to different insights into factors that structure microbial communities.  Applied and Environmental Microbiology 73, 5 (2007), 1576–85.
  * [18] McInnes, L., Healy, J., and Melville, J. UMAP: Uniform manifold approximation and projection for dimension reduction, February 2018.
  * [19] Qu, T., and Cai, Z. A fast multidimensional scaling algorithm.  In 2015 IEEE International Conference on Robotics and Biomimetics (ROBIO), pp. 2569–2574.
  * [20] Reiman, D., Metwally, A. A., Sun, J., and Dai, Y. Popphy-cnn: A phylogenetic tree embedded architecture for convolutional neural networks to predict host phenotype from metagenomic data.  IEEE Journal of Biomedical and Health Informatics 24, 10 (2020), 2993–3001.
  * [21] Tenenbaum, J. B., de Silva, V., and Langford, J. C. A global geometric framework for nonlinear dimensionality reduction.  Science 290, 5500 (2000), 2319–23.
  * [22] van der Maaten, L., and Hinton, G. Visualizing data using t-sne.  Journal of Machine Learning Research 9, 86 (2008), 2579–2605.
  * [23] Witten, D. M., and Tibshirani, R. Supervised multidimensional scaling for visualization, classification, and bipartite ranking.  Computational Statistics & Data Analysis 55, 1 (2011), 789–801.
  * [24] Yang, F., Yang, W., Gao, R., and Liao, Q. Discriminative multidimensional scaling for low-resolution face recognition.  IEEE Signal Processing Letters 25, 3 (2018), 388–392.
  * [25] Yang, T., Liu, J., McMillan, L., and Wang, W. A fast approximation to multidimensional scaling.  In Proceedings of the ECCV Workshop on Computation Intensive Methods for Computer Vision (CIMCV).
  * [26] Zheng, J. X., Pawar, S., and Goodman, D. F. M. Graph drawing by stochastic gradient descent.  IEEE Transactions on Visualization and Computer Graphics 25, 9 (2019), 2738–2748.


# Breaking a Bloch-wave interferometer: quasiparticle species-specific
temperature-dependent nonequilibrium dephasing

Joseph B. Costello  Physics Department, University of California, Santa
Barbara, California 93106, USA  Institute for Terahertz Science and
Technology, University of California, Santa Barbara, California 93106, USA
Seamus D. O’Hara  Physics Department, University of California, Santa Barbara,
California 93106, USA  Institute for Terahertz Science and Technology,
University of California, Santa Barbara, California 93106, USA  Qile Wu
Physics Department, University of California, Santa Barbara, California 93106,
USA  Institute for Terahertz Science and Technology, University of California,
Santa Barbara, California 93106, USA  Moonsuk Jang  Physics Department,
University of California, Santa Barbara, California 93106, USA  Institute for
Terahertz Science and Technology, University of California, Santa Barbara,
California 93106, USA  Loren N. Pfeiffer  Electrical Engineering Department,
Princeton University, Princeton, New Jersey 08544, USA  Ken W. West
Electrical Engineering Department, Princeton University, Princeton, New Jersey
08544, USA  Mark S. Sherwin  Physics Department, University of California,
Santa Barbara, California 93106, USA  Institute for Terahertz Science and
Technology, University of California, Santa Barbara, California 93106, USA

###### Abstract

Recently, high-order sideband polarimetry has been established as an
experimental method that links the polarization of sidebands to an
interference of Bloch wavefunctions. However, the robustness of sideband
polarizations to increasing dephasing remains to be explored. Here, we
investigate the dependence of high-order sideband generation in bulk gallium
arsenide on dephasing by tuning temperature. We find that the intensities of
the sidebands, but not their polarizations, depend strongly on temperature.
Using our polarimetry method, we are able to isolate the contributions of
electron-heavy hole (HH) and electron-light hole (LH) pairs to sideband
intensities, and separately extract the nonequilibrium dephasing coefficients
associated with the longitudinal optical (LO) phonons and acoustic (A) phonons
for each species of electron-hole pair. We find that
$\Gamma_{\text{HH},\text{A}}=6.1\pm 1.6$ $\mu$eV/K,
$\Gamma_{\text{LH},\text{A}}<1.5$ $\mu$eV/K,
$\Gamma_{\text{HH},\text{LO}}=14\pm 3$ meV, and
$\Gamma_{\text{LH},\text{LO}}=30\pm 3$ meV.

††preprint: APS/123-QED

##  I Introduction

One of the chief aims of modern condensed matter physics is to study the
coherent properties of quasiparticles in the form of Bloch waves in
crystalline solids. In contrast with atomic and photonic systems, the dynamics
of coherence in solids is complicated by interactions between the Bloch waves
of interest and other quasiparticles. These interactions can lead to dephasing
which destroys quantum coherence. The importance of dephasing processes in
materials has led to a great deal of interest in understanding the role of
dephasing in crystals in general [1, 2, 3, 4] and in specific systems [5, 6,
7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17].

A powerful optical technique that is deeply affected by dephasing is
interferometry. Interferometry allows experimental access to quantities of
interest such as wavefunction phases and quantum exchange statistics, which
are otherwise difficult to measure [18, 19, 20]. Because interferometric
visibility relies on coherent processes, it can be reduced or extinguished by
dephasing. Despite this difficulty, interference in solids has been observed
in a variety of systems. In some experiments, coherences are protected by
engineering devices that host states with very small dephasing, such as chiral
edge states of two-dimensional electron gases in quantum wells [21, 13, 22,
18]. Other coherent processes are detected over time periods that are
comparable to or shorter than the dephasing time. For example, in high
harmonic generation (HHG) in solids, an intense laser drives charge carriers
and causes the emission of harmonics with frequencies that are multiples of
the driving laser frequency. The driving laser is sufficiently strong that the
harmonics are generated before the carriers are dephased [23, 19, 24, 25, 26].

Interferometric studies of HHG are exciting because HHG is a non-destructive
process that can be applied to a wide variety of solid state systems such as
traditional semiconductors [25], ultra wide gap insulators [27], and
topological insulators [20]. Depending on the materials and driving fields,
different intraband trajectories, multiple interband ionization pathways, or
the presence of both intraband and interband contributions can lead to quantum
interferences [19, 24, 26, 28, 29, 15]. However, the richness of microscopic
interference sources makes it difficult to attribute a given HHG signal to any
particular contribution.

High-order sideband generation (HSG) is a nonlinear optical process similar to
the interband processes in HHG. HSG occurs in semiconductors when a relatively
weak near-infrared (NIR) laser tuned to the bandgap creates an electron-hole
pair, which is then accelerated by a strong low-frequency laser. The
recombination of the electron and hole results in the emission of sideband
photons. Unlike HHG, ballistic recollisions of electron-hole pairs are the
only known coherent processes in HSG [30]. The absence of other mechanisms
makes the analysis of HSG relatively simple. If there are multiple bands that
are degenerate or close in energy, then these bands can all contribute to the
sideband signals. These different contributions interfere, influencing the
polarizations of the sidebands. A newly developed technique called high-order
sideband polarimetry, has been used to reconstruct the Bloch wavefunctions of
holes in bulk gallium arsenide (GaAs) [31], measure the effects of the dynamic
phases of Bloch waves in GaAs [32], and probe Berry curvatures in GaAs quantum
wells [33]. Based on polarimetry of high-order sidebands, methods of directly
extracting the dephasing rates and information about band structures have also
been proposed [34]. HSG has been demonstrated in several different materials
up to room temperature [35, 36, 37, 38]. However, the robustness of sideband
polarization in HSG to increasing dephasing has not been investigated.
Specifically, in bulk GaAs and GaAs-based quantum wells, quantum interferences
in sideband polarizations have been studied only at low temperature, where
dephasing is minimal [30, 39].

Dephasing of interband excitations in GaAs has been studied previously, in
both theory and experiment [40, 41, 42, 43, 44, 45]. Experimentally, this is
normally done by examining the width of optical peaks, for example, in linear
absorption and photoluminescence. These experiments are complicated by the
fact that GaAs has two species of holes, heavy holes (HH) and light holes
(LH). These species of holes are degenerate at the band edge. This makes it
difficult to experimentally distinguish the dephasing rates of their
corresponding electron-hole pairs, and to our knowledge previous experiments
have not been able to separately measure them.

Here, we investigate the role of dephasing in HSG in bulk GaAs by tuning the
temperature from 25 K to 200 K. We perform these experiments in the low-
density regime where interactions between photo-generated electron-hole pairs
can be ignored, so the primary sources of temperature-dependent dephasing come
from interactions with acoustic and longitudinal optical (LO) phonons [46].
Since high-order sideband polarimetry can distinguish between contributions to
sideband polarization from different hole species, we are able to analyze
temperature-dependent sideband polarization data to separately extract the
temperature-dependent dephasing rates of the electron-HH and electron-LH pairs
and determine the dephasing coefficents associated with the acoustic phonons
and LO phonons for each species of electron-hole pair.

Figure 1: Temperature dependent high-order sideband generation (HSG). (a) The
band structure near the $\Gamma$ point ($\textbf{k}=0$) of gallium arsenide
(GaAs). There is an electron (E) band, a heavy-hole (HH) band, and a light-
hole (LH) band. The bandgap is shrunk by several orders of magnitude for
visualization purposes. In HSG, electron-hole pairs are created at the
$\Gamma$ point and accelerated to higher quasimomenta such that when
annihilation occurs sideband photons are emitted at higher energies. (b) Real
space trajectories of electrons and holes that contribute to the same order of
sideband. 1) The E-HH (blue) and E-LH (red) pairs are created by the near-
infrared (NIR) laser, 2) and are accelerated by a strong linearly polarized
terahertz (THz) field. 3) The electrons and holes recombine, releasing a
photon at a higher energy. (c) A sideband spectrum measured at 25 K with the
NIR laser polarized parallel to the THz field. Inset: a cartoon of the
experimental setup. The NIR and THz lasers are focused collinearly on an
epilayer of GaAs mounted onto a sapphire substrate. The THz field is reflected
by an indium tin oxide (ITO) layer which creates an optical enhancement cavity
for the THz field. The NIR laser and sideband fields are transmitted through
the ITO to be measured. (d) Temperature-dependent absorbance spectra from 25 K
(blue) to 200 K (red) in dB, with diamonds representing the NIR photon
energies used in HSG experiments at each temperature. For further details on
this measurement, see Appendix B. Each spectrum is offset by a 1 dB step for
visiblity. e) Sideband intensities as functions of sideband order at different
temperatures ranging from 25 K (blue) to 200 K (red). Experiments were
performed at 25 K ($\hbar\omega_{\text{NIR}}=1.512$ eV), 100 K
($\hbar\omega_{\text{NIR}}=1.500$ eV), 120 K ($\hbar\omega_{\text{NIR}}=1.493$
eV), 140 K ($\hbar\omega_{\text{NIR}}=1.486$ eV), 160 K
($\hbar\omega_{\text{NIR}}=1.477$ eV), 180 K ($\hbar\omega_{\text{NIR}}=1.472$
eV), and 200 K ($\hbar\omega_{\text{NIR}}=1.462$ eV). The intensities are
normalized to that of the 2nd order sideband at 25 K. In (c) and (e), the NIR
laser field is polarized parallel to the THz electric field, and only the
components along the THz polarization are shown. See Appendix A for details.

##  II High-order Sideband Polarimetry

The experiments described here were performed in a bulk gallium arsenide
(GaAs) epilayer. A near-infrared (NIR) laser was used to create electron-hole
pairs and a strong terahertz (THz) laser was used to drive their recollisions.
The THz source used in this paper was the UCSB millimeter-wave Free Electron
Laser (FEL) [47], with a frequency of $449\pm 1$ GHz [48], a field strength of
$65\pm 3$ kV/cm at the GaAs epilayer, and a pulse length of $40$ ns. The
electron-hole pairs were accelerated by the THz field in the (001) plane of
GaAs. The energies of sidebands are

$\hbar\omega_{\text{SB},n}=\hbar\omega_{\text{NIR}}+n\hbar\omega_{\text{THz}},$
(1)

where $\hbar$ is the reduced Planck constant, $n$ is an integer called the
sideband order, $\omega_{\text{SB},n}$ is the frequency of the nth-order
sideband, $\omega_{\text{NIR}}$ is the frequency of the NIR laser, and
$\omega_{\text{THz}}$ is the THz frequency. Because of the inversion symmetry
of the (001) plane in GaAs, $n$ must be even.

Figure 1 (a) shows the relevant band structure of GaAs with one electron (E)
band, a heavy-hole (HH) band, and a light-hole (LH) band [49]. Due to the
degeneracy of the HH and LH bands at the band edge, the NIR laser creates both
E-HH and E-LH pairs, which are then accelerated by the THz field. For a
sideband of sufficiently low energy, there are classical recollision pathways
associated with both the E-HH and E-LH pairs, as shown in Fig. 1 (b).
Depending on when the electron-hole pairs are created, they are annihilated at
different times with different kinetic energies, leading to emission of
different sidebands. Figure 1 (c) shows a typical sideband spectrum, with
sideband orders from $n=2$ to $n=52$ visible. This particular spectrum was
taken at $\lambda_{\text{NIR}}=820$ nm with the NIR laser polarized parallel
to the THz field. The inset of Fig. 1 (c) shows the experimental setup of our
experiments, with the NIR and THz lasers collinearly focused on a $500$ nm
epitaxially grown GaAs that was mounted onto a sapphire substrate. An indium
tin oxide (ITO) film, which reflects the THz field to create an enhancement
cavity and is transmissive to the NIR laser and sideband fields, was deposited
onto the opposite side of the substrate. The sample is the same one as used in
Ref. [32] and is prepared in the same way as described in Ref. [31]. On top of
the ITO film, there is a silicon dioxide layer acting as an anti-reflection
coating for the NIR laser and sideband fields.

Each sideband contains contributions from the two species of holes, which have
different angular momenta. Since the total angular momentum must be conserved
in the electron-photon interaction, this difference results in different
photon helicity produced by different types of electron-hole pair. Therefore,
sideband polarizations can be thought of as resulting from interferences
between the different quantum paths followed by the E-LH and E-HH pairs.

The falloff of sideband intensity with increasing sideband order is controlled
by dephasing. Previous work has established that HSG in bulk GaAs persists up
to 170 K [35], but no previous study has investigated the polarizations of
these sidebands at temperatures above 60K. The analysis of HSG is greatly
simplified by tuning the NIR laser to the bandgap $E_{\rm g}$, i.e.,
$\hbar\omega_{\text{NIR}}=E_{\rm g}$ [31]. In this way all electrons and holes
are created at quasimomentum ${\bf k}=0$. However, the bandgap decreases with
increasing temperature. Figure 1 (d) shows the optical absorption of our
sample at different temperatures. The peaks get broader with increasing
temperature, which is indicative of stronger near-equilibrium dephasing at
higher temperatures.

We performed high-order sideband polarimetry at the temperatures listed in
Fig. 1 (d) and Fig. 1(e). In HSG experiments, the NIR laser was tuned to the
temperature-dependent bandgap of GaAs at each temperature. Figure 1 (e) shows
the sideband intensity as a function of sideband order for all temperatures we
investigated. For higher temperatures, one expects that the stronger dephasing
will result in sideband intensities that are generally lower, and hence fewer
sidebands will be observed. However, below 160 K, we observed similar numbers
of sidebands. This is encouraging for applications of HSG in broadband
frequency combs [48]. Our results indicate that cooling down to below liquid-
nitrogen temperatures is not necessary to obtain acceptable sideband signals.

##  III Temperature Dependence of Sideband Polarization

Figure 2: Temperature-dependent HSG polarimetry. The linear orientation angle
$\alpha$ (upper panel) and ellipticity angle $\gamma$ (lower panel) of
sidebands are plotted as functions of sideband order for temperatures ranging
from 25 K (blue) to 200 K (red). The colors match those in Fig. 1 (d) and (e).
The angles $\alpha$ and $\gamma$ are defined with respect to the THz field
(cartoon on the right). For each temperature, data for four linear
polarizations of the NIR laser ($\gamma_{\text{NIR}}=0^{\circ}$) are shown,
with $\alpha_{\text{NIR}}=90^{\circ}$ (diamonds),
$\alpha_{\text{NIR}}=45^{\circ}$ (triangles), $\alpha_{\text{NIR}}=0^{\circ}$
(circles), and $\alpha_{\text{NIR}}=-45^{\circ}$ (squares). Above 140 K, fewer
sidebands are measured. At 200 K, no sidebands were measured for
$\alpha_{\text{NIR}}=90^{\circ}$. See Appendix A for details on polarization
measurement.

We measure the polarizations of sidebands using the method of Stokes
polarimetry (see Appendix A for details). This method is sensitive to both the
linear orientation angle $\alpha$ and the ellipticity angle $\gamma$ (see the
polarization ellipse on the right-hand side of Fig. 2). Numerically, we
constrain these angles to be in the ranges $-90^{\circ}\leq\alpha\leq
90^{\circ}$ and $-45^{\circ}\leq\gamma\leq 45^{\circ}$. A linear polarization
corresponds to $\gamma=0$, and right-handed (left-handed) circularly polarized
light corresponds to $\gamma=45^{\circ}$ ($\gamma=-45^{\circ}$).

Figure 2 shows the polarimetry results for all temperatures we investigated.
At each temperature, four polarimetry scans were performed for four different
linear polarizations of the NIR laser. The four linear orientation angles of
the NIR laser were $\alpha_{\text{NIR}}=90^{\circ}$ (diamonds),
$\alpha_{\text{NIR}}=45^{\circ}$ (triangles), $\alpha_{\text{NIR}}=0^{\circ}$
(circles), and $\alpha_{\text{NIR}}=-45^{\circ}$ (squares). The upper and
lower panels show the $\alpha$ and $\gamma$ angles for the sidebands. The
sideband polarizations exhibit remarkable robustness to temperature, with
small variations up to 160 K. Above 160 K, reduced signal-to-noise ratio leads
to more deviation in higher-order sidebands, but even at 200 K the
polarizations of the lower-order sidebands are very similar to those at lower
temperatures.

These results reinforce our understanding that the sideband polarizations are
set by the interference of Bloch waves. In our HSG experments, this
interference arises from different dynamic phases accumulated by the E-HH
(E-LH) pair from the creation time $t_{0}$ to annihilation time $t_{f}$,

$A_{\text{HH}(\text{LH})}(t_{f},t_{0})=-\frac{1}{\hbar}\int_{t_{0}}^{t_{f}}dt^{\prime}E_{\text{HH}(\text{LH})}(t^{\prime}).$
(2)

Here, $t_{0}$ and $t_{f}$ depend on the sideband order as well as the species
of electron-hole pairs. The relative energy of the E-HH and E-LH pairs,
$E_{\rm HH}$ and $E_{\rm LH}$, contain the bandgap $E_{\rm g}$, which is the
only temperature-dependent factor in the dynamic phases. Since the NIR laser
is tuned to the bandgap at each temperature, this temperature dependence is
not relevant for the sideband polarization [31]. The dephasing, however, does
depend strongly on temperature. This suggests that the electron-hole
coherences are damped faster at higher temperature by the stronger dephasing,
resulting in the lower sideband intensities, as shown in Fig. 1 (e). In
addition, the electron-hole pairs that survive until annihilation and sideband
emission experience the same coherent acceleration, and therefore the
information contained in the sideband polarization is nearly independent of
temperature.

##  IV Extracting Dephasing Rates from Sideband Polarimetry

Figure 3: Temperature-dependent falloffs of the E-HH and E-LH contributions to
sideband intensity. (a) The absolute values of the propagators for the E-HH
and E-LH pairs, $\varsigma_{\text{HH}}$ (left, squares) and
$\varsigma_{\text{LH}}$ (right, diamonds) for different temperatures running
from 25 K (blue) to 200 K (red) as functions of sideband order. The colors
match those in Fig. 2. (b) The absolute values of $\varsigma_{\text{HH}}$
(left, circles) and $\varsigma_{\text{LH}}$ (right, triangles) of different
sideband orders as functions of temperature. Each plot is normalized to the 25
K data to isolate the dependence of the falloff for each sideband on
temperature. The dotted lines in each frame show the results from fitting all
data for each hole species to Eq. 6. For each fit, there are two free
parameters, $\Gamma_{j,\text{A}}$ and $\Gamma_{j,\text{LO}}$, where $j$ labels
the hole species. Different sideband orders are associated with different
acceleration time $\tau_{j}(n)$ for the recollision processes, so the dotted
lines for different sideband orders are slightly different despite the
dephasing function $\Gamma_{j}(T)$ being the same.

Although the polarization of sidebands does not depend strongly on
temperature, the intensities of these sidebands can carry information on the
temperature dependence of dephasing. By using the results of Ref. [31], the
temperature dependence of the dephasing rates of the E-LH and E-HH pairs can
be extracted separately. The electric field of the $n$th-order sideband can be
written as

$\mathbb{P}_{n}\propto\sum_{s}\begin{pmatrix}{\bf D}_{\text{E-HH},s}\\\ {\bf
D}_{\text{E-LH},s}\end{pmatrix}^{\dagger}\begin{pmatrix}\varsigma_{\text{HH},n}&0\\\
0&\varsigma_{\text{LH},n}\end{pmatrix}\begin{pmatrix}{\bf
D}_{\text{E-HH},s}\\\ {\bf D}_{\text{E-LH},s}\end{pmatrix}\cdot{\bf
F}_{\text{NIR}}$ (3)

where $s$ labels the two-fold degeneracy in the electron-hole states, ${\bf
D}_{\text{E-HH (E-LH)},s}$ is the electric dipole vector associated with the
E-HH (E-LH) pair labeled by $s$, ${\bf F}_{\text{NIR}}$ is the electric field
of the NIR laser, and $\varsigma_{\text{HH}(\text{LH}),n}$ is the propagator
of the E-HH (E-LH) pair. The dynamics of electron-hole recollisions is
captured by the propagator $\varsigma_{j,n}$ ($j={\rm HH,LH}$), which assumes
the form,

$\varsigma_{j,n}\propto\int_{0}^{2\pi/\omega_{\rm THz}}dte^{in\omega_{\rm
THz}t}\int_{-\infty}^{t}dt^{\prime}e^{iA_{j}(t,t^{\prime})-(t-t^{\prime})\Gamma_{j}(T)/\hbar},$
(4)

where $A_{j}(t,t^{\prime})$ is the dynamic phase defined by Eq. 2, and
$\Gamma_{j}(T)/\hbar$ is the dephasing rate written as a function of
temperature $T$. The first integral is a Fourier transform to isolate the
electron-hole pairs which have gained kinetic energy $n\hbar\omega_{\rm THz}$
upon recombination and therefore contribute to the $n$th order sideband. The
second integral sums over the contributions from electron-hole pairs created
at different times $t^{\prime}$. In this expression, only the recollision
pathways associated with electron-hole pairs created at zero quasimomentum are
included. Because the NIR laser in our HSG experiment is relatively weak, the
population of electrons and holes is not high enough for carrier-carrier
scattering effects to be an important source of dephasing. We ensure that this
is the case by confirming that the sideband intensities are linearly
proportional to the NIR laser power (we used 100 mW in these experiments). We
model the dephasing function $\Gamma_{j}(T)$ ($j={\rm HH,LH}$) as

$\Gamma_{j}(T)=\Gamma_{0}+\Gamma_{j,\text{A}}T+\Gamma_{j,\text{LO}}n_{\text{LO}}(T),$
(5)

which includes a temperature-independent dephasing constant $\Gamma_{0}$, and
contributions from the long-wavelength acoustic phonons and the longitudinal
optical (LO) phonons proportional to the phonon occupation numbers. We assume
that the occupation number for the long-wavelength acoustic phonons is
proportional to temperature, while the occupation number for the long-
wavelength LO phonons is taken as Bose-Eistein distribution, $n_{\rm
LO}(T)=1/[\exp(\hbar\omega_{\rm LO}/k_{\rm B}T)-1]$, where $\hbar\omega_{\rm
LO}=36.6$ meV is the energy of the LO phonons [39] and $k_{\rm B}$ is the
Boltzmann constant. We refer to $\Gamma_{j,\text{A}}$ and
$\Gamma_{j,\text{LO}}$ as dephasing coefficients since both need to be
multiplied by additional factors to give the dephasing rate.

Figure 3 (a) shows the absolute values of the propagators,
$|\varsigma_{\text{HH},n}|$ and $|\varsigma_{\text{LH},n}|$, as functions of
sideband order. As expected, the propagators have larger amplitudes at lower
temperatures. Intriguingly, the falloffs of the propagators for the E-HH and
E-LH pairs are slightly different. This is due to slight differences in
dephasing rates, creation times, and annihilation times for the E-HH and E-LH
pairs.

We assume here that the dephasing rate $\Gamma_{j}/\hbar$ is independent of
quasimomentum, and can be understood as the mean dephasing rate experienced by
the electron-hole pairs during acceleration from the THz field. This
assumption implies that the dephasing rate is constant in time for a given
electron-hole pair, which allows us to isolate the dephasing dependence in
$|\varsigma_{j,n}|$ as

$|\varsigma_{j,n}|\propto e^{-\Gamma_{j}(T)\tau_{j}(n)/\hbar},$ (6)

where $\tau_{j}$ is the time between the creation and annihilation of the
electron-hole pair labeled by $j$ in the shortest classical recollision
pathway (referred to as the acceleration time of the electron-hole pairs). We
have assumed that the dephasing investigated here is sufficiently strong such
that, for each species of the electron-hole pair, there is one shortest
recollision pathway that dominantly determines the polarization of a sideband.
To simplify the analysis, we consider here classical recollisions, in which
the electrons and holes are created and annihilated at the same position. In
our HSG experiment, the average kinetic energy gain of the electron-hole pairs
in a THz cycle, namely, the ponderomotive energy $U_{{\rm
p},j}={e^{2}F^{2}_{\rm THz}}/{4\mu_{j}\omega_{\rm THz}^{2}}$ ($j={\rm
HH,LH}$), is much larger than the sideband offset energies and the dephasing
function $\Gamma_{j}(T)$. Here, $e$ is the elementary charge, $F_{\rm THz}$ is
the THz field strength, and $\mu_{\rm
HH}=[m_{c}^{-1}+\gamma_{1}-2\gamma_{2}|{\bf n}|]^{-1}$ ($\mu_{\rm
LH}=[m_{c}^{-1}+\gamma_{1}+2\gamma_{2}|{\bf n}|]^{-1}$) is the reduced mass of
the E-HH (E-LH) pairs along the linear polarization of the THz field. The
parameter $m_{c}=0.067m_{0}$ is the effective mass of the conduction band
($m_{0}$ is the electron rest mass), $\gamma_{1}=6.98$, $\gamma_{2}=2.06$ are
two Luttinger parameters [50], and the vector ${\bf
n}=\left(({\sqrt{3}}/2)\sin 2\theta,-({\sqrt{3}\gamma_{3}}/{2\gamma_{2})}\cos
2\theta,-{1}/{2}\right)$ is defined by the Luttinger parameters $\gamma_{2}$
and $\gamma_{3}=2.93$ [50], together with the angle between the linear
polarization of the THz field and the [110] crystal direction of GaAs,
$\theta=67^{\circ}$. As discussed in Ref. [34], in this case, the THz field
can be approximated as linear in time during the recollision processes of the
electron-hole pairs. With this approximation, the acceleration time associated
with the $n$th-order sideband can be explicitly calculated as [34, 32]

$\displaystyle\tau_{j}(n)=\frac{\sqrt{3}}{\omega_{\rm
THz}}\left(\frac{2n\hbar\omega_{\rm THz}}{U_{{\rm p},j}}\right)^{1/4},$ (7)

The acceleration times $\tau_{j}$ are different for the E-HH and E-LH pairs
because of their different reduced masses, leading to different trajectories
that cause quantum interferences.

Figure 3 (b) shows $|\varsigma_{j,n}|$ as a function of temperature for
sideband orders from $n=10$ to $n=20$. For each sideband order, the values of
$|\varsigma_{j,n}|$ are normalized to its value at 25 K. This normalization
allows all of the data to be fit by exponential functions (Eq. 6) with only 2
free parameters for each hole species, $\Gamma_{j,\text{A}}$ and
$\Gamma_{j,\text{LO}}$. The dotted lines in Fig. 3 (b) show the fitting
results. There are different lines for different sideband orders because of
the difference in the acceleration time $\tau_{j}(n)$.

Figure 4: Extracted coefficients for the temperature dependent part of
dephasing rates. The coefficients $\Gamma_{j,\text{A}}$ and
$\Gamma_{j,\text{LO}}$ in the dephasing function (Eq. 5) are extracted by
fitting to scaled $|\varsigma_{\text{HH}}|$ and $|\varsigma_{\text{LH}}|$
shown in Fig. 3 (b) with Eq. 6. The red (blue) solid horizontal lines
represent the results from fitting to the data for all sideband orders with a
dephasing function $\Gamma_{\text{LH}}(T)$($\Gamma_{\text{HH}}(T)$) that is
assumed to be independent of the sideband order. The numerical values are
$\Gamma_{\text{HH},\text{A}}=6.1\pm 1.6$ $\mu$eV/K,
$\Gamma_{\text{LH},\text{A}}<1.5$ $\mu$eV/K,
$\Gamma_{\text{HH},\text{LO}}=14\pm 3$ meV, and
$\Gamma_{\text{LH},\text{LO}}=30\pm 3$ meV, which produce the dotted lines in
Fig. 3 (b). These error bars are the standard deviations by taking the square
roots of the variance of the fits.

Figure 4 shows the parameters $\Gamma_{j,\text{A}}$ and $\Gamma_{j,\text{LO}}$
that are returned by the fitting procedure. The diamonds (dots) refer to fits
for each sideband order individually for the E-HH (E-LH) pairs, and the solid
lines show the results from fitting to the data for all sideband orders
simultaneously. The results of the fit are $\Gamma_{\text{HH},\text{A}}=6.1\pm
1.6$ $\mu$eV/K, $\Gamma_{\text{LH},\text{A}}<1.5$ $\mu$eV/K,
$\Gamma_{\text{HH},\text{LO}}=14\pm 3$ meV, and
$\Gamma_{\text{LH},\text{LO}}=30\pm 3$ meV. The errors reported are the square
roots of the variances of the fit. The fits are likely less reliable for the
dephasing coefficents associated with the acoustic phonons since the
contribution from the LO phonon dominates for most of the sampled temperature
range. This result agrees with previous results that the LO phonon is the
dominant source of dephasing in bulk GaAs above 80 K [42, 39].

It is important to note that the coefficients $\Gamma_{j,\text{LO}}$ must be
multiplied by the LO phonon occupation number to give their contribution to
$\Gamma_{j}(T)$ in Eq. 5. Thus, for example, at 100 K, where the phonon
occupation number is 0.015, the LO-phonon contributions to $\Gamma_{j}(T)$ are
$\Gamma_{\rm HH,LO}n_{\rm LO}=0.21$ meV and $\Gamma_{\rm HH,LO}n_{\rm
LO}=0.45$ meV. For comparison, $\Gamma_{j}(T=35{\rm K})$, which is close to
$\Gamma_{0}$ because the LO phonon occupation is near 0, has been estimated to
be $4.8\hbar\omega_{\rm THz}=9$ meV [32].

##  V Discussion

Our results are not the first to measure dephasing in GaAs, but our method
based on sideband polarimetry represents a new way to measure dephasing that
differs from previous experiments. The ability of our method to distinguish
the contributions from different Bloch waves is a key advantage. Although
theory has long predicted that the scattering and dephasing rates should
differ for different hole species [46, 43], this is to our knowledge the first
experiment that differentiates the dephasing coefficients of these holes near
the bandgap. Previous investigations of dephasing in GaAs have used the
linewidth of photoluminescence or absorption spectra, and have been
complicated by the degeneracy of the LH and HH band at the band edge [43, 42].
These experiments found $\Gamma_{\text{A}}=13\pm 3$ $\mu$eV/K and
$\Gamma_{\text{LO}}=30.4\pm 4$ meV, but they were unable to differentiate HH
from LH contributions, and had little confidence in the value of
$\Gamma_{\text{A}}$ [42]. These methods have particular difficulty at high
temperatures, where strain-induced splittings of the HH and LH bands are
washed out by broadened peaks. Four-wave mixing has also been used to probe
dephasing in GaAs [45], but is similarly only able to examine low-temperature
physics. It is interesting that these near-equilibrium measurements seem to
agree with our measurement of $\Gamma_{\text{LH},\text{LO}}$. This is perhaps
due to the fact that the LH peak is broader due to higher dephasing, and could
dominate these linewidth measurements.

Theoretical explorations of $\Gamma_{\text{LO}}$ in GaAs have found that
dephasing of excitons depends heavily on the effective mass of the hole, with
lower hole masses corresponding to higher dephasing rates [43]. This study
found $\Gamma_{\text{HH},\text{LO}}=21$ meV for excitons in near equilibrium.
The difference from our results could come from the fact that the dephasing
experienced by the electron-hole pairs in our HSG experiment are in continuum
states far from equilibrium. Indeed, the methods in Ref. [43] provide insights
only into the near-equilibrium dephasing. Recently, dephasing in systems far
from equilibrium in the context of HHG has been studied theoretically [16],
including effects of phonons [17]. However, to our knowledge, species-specific
nonequilibrium dephasing rates in semiconductors have not been previously
measured. Our results on GaAs pave the way to measuring the nonequilibrium
dephasing rates of a large class of semiconducting materials.

A further advantage of our method is the possibility of measuring the
dependence of dephasing on the direction of electron-hole acceleration. By
rotating the THz polarization with respect to the crystal, the electron-hole
pairs will be driven along different directions. By studying this angular
dependence of sideband polarization, it is possible to probe dephasing in
different crystal directions. This sensitivity is not possible in linear
optical methods such as linear absorption. Future experiments could be
particularly useful in probing the directional dephasing of systems that have
significant anisotropy. In GaAs, changing the direction of electron-hole
acceleration could probe the dependence of dephasing on effective masses of
holes, since the effective masses of both the HH and LH depend on crystal
orientation.

##  VI Conclusion

As application of the quantum effects in materials advances, experimental
knowledge of the coherent properties of Bloch wavefunctions becomes more
critical. Polarimetry of high-order sidebands has been shown as an exciting
interferometry method that provides insight into the phases of Bloch
wavefunctions and allows extraction of key quantities in electronic structures
[31, 32, 34, 38, 37, 36]. The fact that these polarization states are robust
over a wide range of temperatures is encouraging both to the viability of
further low cost HSG experiments in GaAs and to the extension of HSG to other
materials which have greater dephasing rates. Furthermore, in devices that
take advantage of coherent quantum properties, knowledge of dephasing is
critical. We have demonstrated that temperature-dependent high-order sideband
polarimetry can be used to differentiate the temperature-dependent dephasing
of different species of electron-hole pairs, even if they are very close in
energy. Future experiments can complete this analysis by examining the
temperature-independent contribution to dephasing, and by extending this
method to other materials that host more exotic wavefunctions.

###### Acknowledgements.

We acknowledge N. Agladze for assistance in maintaining and operating the UCSB
millimeter-wave free electron laser. This work was funded by NSF DMR 2004995.
J.B. Costello and M.S. Sherwin conceptualized this work. Q. Wu produced the
theory. J.B. Costello, S.D. O’hara, and M. Jang performed the experiments.
J.B. Costello performed data analysis. L.N. Pfeiffer and K.W. West grew the
GaAs crystal. M.S. Sherwin acquired funding and administered the project.

##  Appendix A  Stokes Polarimetry

To measure sideband polarizations, we used a method called Stokes polarimetry,
following the same procedure as described in Ref. [31]. Stokes polarimetry was
performed by measuring the intensity of a light wave after it propagated
through a rotating quarter-wave plate and a linear polarizer as a function of
the rotation angle of the quarter-wave plate. In our experiments, the
transmission axis of the linear polarizer was parallel to the linear
polarization of the THz field defined as horizontal. The intensity $I$ of the
outgoing light can be written as

$\displaystyle I(\varphi)=$
$\displaystyle\frac{S_{0}}{2}+\frac{S_{1}}{4}[1+\cos(4\varphi)]$
$\displaystyle+\frac{S_{2}}{4}\sin(4\varphi)-\frac{S_{3}}{2}\sin(2\varphi),$
(8)

where $\varphi$ is the angle between the fast axis of the rotating quarter-
wave plate and the horizontal, and $S_{0}$, $S_{1}$, $S_{2}$, and $S_{3}$ are
the Stokes parameters, which define the polarization state of the light. We
sampled 16 different rotation angles of the quarter-wave plate to extract the
Stokes parameters of the sidebands. The Stokes parameters are related to the
linear orientation angle $\alpha$ and ellipticity angle $\gamma$ through the
following equations,

$\tan(2\alpha)=\frac{S_{2}}{S_{1}},$ (9)
$\tan(2\gamma)=\frac{S_{3}}{\sqrt{S_{1}^{2}+S_{2}^{2}}}.$ (10)

Because our experiments always use the Stokes polarimeter, the full
intensities of the sidebands described by the Stokes parameter $S_{0}$ were
not directly recorded. The spectra displayed in Fig. 1 (e) are actually the
intensities $I(\varphi)$ of the sidebands at $\varphi=0^{\circ}$ for a
horizontal polarization of the NIR laser ($\alpha_{\text{NIR}}=0^{\circ}$).
They should be understood as the intensities of the horizontal components of
the sidebands. The intensity $I(\varphi)$ of a sideband tends to be the
strongest when the NIR laser is horizontally polarized, because in this case
the sidebands tend to be close to horizontally polarized and have maximum
transmission through the linear polarizer in the Stokes polarimeter.
Similarly, for a vertical polarization of the NIR laser
($\alpha_{\text{NIR}}=90^{\circ}$), the sidebands tend to have minimal
transmission through the linear polarizer, which lowers the signal-to-noise
ratio. For this reason, no sidebands were detected for the vertical
polarization of the NIR laser ($\alpha_{\text{NIR}}=90^{\circ}$) at 200 K.

##  Appendix B  Wavelength Calibration

The sideband intensity data presented here were recorded with two different
detectors, a Hamamatsu R7400U-20 photomultiplier tube (PMT) and an Andor
NewtonEM electron-multiplying charge coupled device (EMCCD). The PMT provides
a time-resolved, single-pixel measurement of a narrow frequency range, while
the EMCCD provides a simultaneous time-integrated measurement in many pixels,
enabling an entire HSG spectrum to be recorded for each NIR laser pulse. The
NIR laser was always much stronger than the sidebands, so the EMCCD tends to
be saturated when we tried to measure the lower-order sidebands, which were
very close to the NIR laser line. For this reason, we use the PMT to measure
the lower-order sidebands and the EMCCD to measure higher-order sidebands. The
Stokes polarimeter was only used with the EMCCD, so the polarizations for the
low-order sidebands ($n<10$) were not measured in this work. See Ref. [31] for
more details.

The quantum efficiencies of both the PMT and EMCCD vary with wavelength over
the range used in this paper. Therefore, to compare HSG spectra at different
temperatures, one needs to calibrate the quantum efficiencies of these two
detectors to account for this wavelength dependence. The wavelengths and
photon energies of the NIR laser we used are listed in Table 1,

Table 1:  The wavelengths and photon energies of the NIR laser used in the HSG experiments at each temperature.  Temperature (K) | Wavelength (nm) | Photon energy (eV)
---|---|---
25 | 820.0 | 1.512
100 | 826.8 | 1.500
120 | 830.2 | 1.493
140 | 834.2 | 1.486
160 | 839.6 | 1.477
180 | 842.5 | 1.472
200 | 848.0 | 1.462

The wavelength of the NIR laser was chosen to be slightly redshifted from the
absorbance peaks in Fig. 1 (d) for a given temperature in order to ensure that
the electron-hole pairs were created at the band edge. The absorbance spectra
shown in Fig. 1 (d) were determined by measuring a broad spectrum of a white
light source with and without the GaAs sample in the transmission path. The
absorbance is calculated (in dB) as

$\text{Abs}=-10~{}\text{log}_{10}(I_{S}/I_{0})$ (11)

where $I_{S}$ ($I_{0}$) is the measured intensity with (without) the GaAs
sample inserted into the optics.

Sidebands up to order 30 were used in polarization measurements, so the
wavelengths measured in this work (including the sidebands and NIR laser)
range from 790 nm to 848 nm. To calibrate both detectors, we used a halogen
white light source and a small tunable monochromator to select a wavelength
range in five nm increments from 778 nm to 848 nm. The monochromator produced
spectra with a wavelength range of roughly 30 nm. The power of the outgoing
light from the monochromator was measured by using a frequency-independent
FieldMaster power meter from Coherent lasers and compared with the signal of
the PMT or EMCCD to determine the detector sensitivity at the selected
wavelength. Since we used the same power of the NIR laser in all sideband
measurements, knowing the detector sensitivities allows for direct comparison
of sideband intensities measured at different temperatures.

The quantum efficiencies of both the PMT and EMCCD decrease with increasing
wavelength, which leads to lower sensitivites at higher temperatures. This
explains why our noise floor is higher for higher temperatures in Fig. 1 (e).

##  Appendix C  Data Analysis

We performed high-order sideband polarimetry for four different linear
polarization of the NIR laser with $\alpha_{\text{NIR}}=90^{\circ}$,
$\alpha_{\text{NIR}}=45^{\circ}$, $\alpha_{\text{NIR}}=0^{\circ}$, and
$\alpha_{\text{NIR}}=-45^{\circ}$, respectively. These measurements allow us
to determine each of the dynamical Jones Matrices up to a constant factor
[33]. The dynamical Jones matrices link the electric fields of the sidebands,
to the electric field of the NIR laser through the following equation,

$\begin{pmatrix}E_{+,n}\\\
E_{-,n}\end{pmatrix}=\begin{pmatrix}T_{++,n}&T_{+-,n}\\\
T_{-+,n}&T_{--,n}\end{pmatrix}\begin{pmatrix}E_{+,\rm NIR}\\\ E_{-,\rm
NIR}\end{pmatrix},$ (12)

where $E_{\pm,n}$ and $E_{\pm,\rm NIR}$ are the components with helicities
$\pm 1$ for the $n$th-order sideband and the NIR laser, respectively, and
$T_{\pm\pm,n}$ are the dynamical Jones matrix elements. Following Ref. [31],
from Eq. 3, the dynamical Jones matrix elements can be written as

$\displaystyle T_{++,n}=T_{--,n}$ $\displaystyle=\frac{2}{3}(\varsigma_{{\rm
HH},n}+\varsigma_{{\rm LH},n})+\frac{n_{z}}{3}(\varsigma_{{\rm
HH},n}-\varsigma_{{\rm LH},n}),$ (13) $\displaystyle T_{+-,n}$
$\displaystyle=\frac{n_{x}+in_{y}}{\sqrt{3}}(\varsigma_{{\rm
HH},n}-\varsigma_{{\rm LH},n}),$ (14) $\displaystyle T_{-+,n}$
$\displaystyle=\frac{n_{x}-in_{y}}{\sqrt{3}}(\varsigma_{{\rm
HH},n}-\varsigma_{{\rm LH},n}).$ (15)

Here, $(n_{x},n_{y},n_{z})$ is a unit vector of the vector

${\bf n}=\left(\frac{\sqrt{3}}{2}\sin
2\theta,-\frac{\sqrt{3}\gamma_{3}}{2\gamma_{2}}\cos
2\theta,-\frac{1}{2}\right),$ (16)

which is determined by the Luttinger parameters $\gamma_{3}=2.93$ and
$\gamma_{2}=2.06$ [50], and the angle between the linear polarization of the
THz field and the [110] crystal direction of GaAs, $\theta=67^{\circ}$. Given
the dynamical Jones matrix elements, we can then calculate the propagators
$\varsigma_{{\rm HH},n}$ and $\varsigma_{{\rm LH},n}$ as

$\displaystyle\varsigma_{\text{LH},n}=\frac{3}{4}\left(T_{++,n}-\frac{2+n_{z}}{\sqrt{3}(n_{x}+in_{y})}T_{+-,n}\right)$
(17)
$\displaystyle\varsigma_{\text{HH},n}=\frac{3}{4}\left(T_{++,n}+\frac{2-n_{z}}{\sqrt{3}(n_{x}+in_{y})}T_{+-,n}\right).$
(18)

To establish confidence intervals in our polarization measurements, we
measured the intensity $I(\varphi)$ of a given sideband 4 times for each
rotation angle $\varphi$ of the quarter-wave plate in the Stokes polarimetry.
The standard deviation calculated from these intensities gives an uncertainty
for each sideband intensity $I(\varphi)$. The uncertainties in the intensity
$I(\varphi)$ propagated throughout the analysis to give the uncertainties in
the angles $\alpha$ and $\gamma$. To determine standard deviations of the
dynamical Jones matrix elements, we use a Monte Carlo method to randomly
sample the angles $\alpha$ and $\gamma$ 10,000 times with their variances set
by the Stokes polarimetry [31]. These deviations then propagate through Eqs.
17 and 18 to determine the uncertainties in the propagators,
$\varsigma_{\text{LH}}$ and $\varsigma_{\text{HH}}$.

## References

  * Skinner [1988] J. L. Skinner, Annual Review of Physical Chemistry 39, 463 (1988).
  * Hsu and Skinner [1984] D. Hsu and J. L. Skinner, The Journal of Chemical Physics 81, 5471 (1984).
  * Hsu and Skinner [1985a] D. Hsu and J. L. Skinner, The Journal of Chemical Physics 83, 2107 (1985a).
  * Hsu and Skinner [1985b] D. Hsu and J. L. Skinner, The Journal of Chemical Physics 83, 2097 (1985b).
  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Physical Review Letters 95, 046404 (2005).
  * Gutman  _et al._ [2008] D. B. Gutman, Y. Gefen, and A. D. Mirlin, Physical Review Letters 101, 126802 (2008).
  * Bagrets  _et al._ [2008] D. A. Bagrets, I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Semiconductors 42, 994 (2008).
  * Konrad  _et al._ [2013] A. Konrad, F. Wackenhut, M. Hussels, A. J. Meixner, and M. Brecht, The Journal of Physical Chemistry C 117, 21476 (2013).
  * Sekine  _et al._ [2003] N. Sekine, Y. Shimada, and K. Hirakawa, Applied Physics Letters 83, 4794 (2003).
  * Gammon  _et al._ [1995] D. Gammon, S. Rudin, T. L. Reinecke, D. S. Katzer, and C. S. Kyono, Physical Review B 51, 16785 (1995).
  * Zhang  _et al._ [1995] Y. Zhang, M. D. Sturge, K. Kash, B. P. van der Gaag, A. S. Gozdz, L. T. Florez, and J. P. Harbison, Superlattices and Microstructures 17, 201 (1995).
  * Rudin and Reinecke [1990] S. Rudin and T. L. Reinecke, Physical Review B 41, 3017 (1990).
  * Neder  _et al._ [2007a] I. Neder, F. Marquardt, M. Heiblum, D. Mahalu, and V. Umansky, Nature Physics 3, 534 (2007a).
  * Yamauchi  _et al._ [2009] Y. Yamauchi, M. Hashisaka, S. Nakamura, K. Chida, S. Kasai, T. Ono, R. Leturcq, K. Ensslin, D. C. Driscoll, A. C. Gossard, and K. Kobayashi, Journal of Physics: Conference Series 193, 012045 (2009).
  * Vampa  _et al._ [2014] G. Vampa, C. McDonald, G. Orlando, D. Klug, P. Corkum, and T. Brabec, Physical Review Letters 113, 073901 (2014).
  * Wang and Du [2021] G. Wang and T.-Y. Du, Physical Review A 103, 063109 (2021).
  * Du and Ma [2022] T.-Y. Du and C. Ma, Physical Review A 105, 053125 (2022).
  * Neder  _et al._ [2007b] I. Neder, N. Ofek, Y. Chung, M. Heiblum, D. Mahalu, and V. Umansky, Nature 448, 333 (2007b).
  * Du  _et al._ [2018] T.-Y. Du, D. Tang, and X.-B. Bian, Physical Review A 98, 063416 (2018).
  * Schmid  _et al._ [2021] C. P. Schmid, L. Weigl, P. Grössing, V. Junk, C. Gorini, S. Schlauderer, S. Ito, M. Meierhofer, N. Hofmann, D. Afanasiev, J. Crewse, K. A. Kokh, O. E. Tereshchenko, J. Güdde, F. Evers, J. Wilhelm, K. Richter, U. Höfer, and R. Huber, Nature 593, 385 (2021).
  * Neder and Marquardt [2007] I. Neder and F. Marquardt, New Journal of Physics 9, 112 (2007).
  * Neder  _et al._ [2007c] I. Neder, M. Heiblum, D. Mahalu, and V. Umansky, Physical Review Letters 98, 036803 (2007c).
  * Jin  _et al._ [2018] J.-Z. Jin, H. Liang, X.-R. Xiao, M.-X. Wang, S.-G. Chen, X.-Y. Wu, Q. Gong, and L.-Y. Peng, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 16LT01 (2018).
  * Hohenleutner  _et al._ [2015] M. Hohenleutner, F. Langer, O. Schubert, M. Knorr, U. Huttner, S. W. Koch, M. Kira, and R. Huber, Nature 523, 572 (2015).
  * Schubert  _et al._ [2014] O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, Nature Photonics 8, 119 (2014).
  * Kim  _et al._ [2019] Y. W. Kim, T.-J. Shao, H. Kim, S. Han, S. Kim, M. Ciappina, X.-B. Bian, and S.-W. Kim, ACS Photonics 6, 851 (2019).
  * Luu  _et al._ [2015] T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T. Hassan, and E. Goulielmakis, Nature 521, 498 (2015).
  * Vampa  _et al._ [2015] G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, and P. B. Corkum, Nature 522, 462 (2015).
  * Ndabashimiye  _et al._ [2016] G. Ndabashimiye, S. Ghimire, M. Wu, D. A. Browne, K. J. Schafer, M. B. Gaarde, and D. A. Reis, Nature 534, 520 (2016).
  * Zaks  _et al._ [2012] B. Zaks, R. B. Liu, and M. S. Sherwin, Nature 483, 580 (2012).
  * Costello  _et al._ [2021] J. B. Costello, S. D. O’Hara, Q. Wu, D. C. Valovcin, L. N. Pfeiffer, K. W. West, and M. S. Sherwin, Nature 599, 57 (2021).
  * O’Hara  _et al._ [2023] S. D. O’Hara, J. B. Costello, Q. Wu, K. West, L. Pfeiffer, and M. S. Sherwin, Bloch-Wave Interferometry of Driven Quasiparticles in Bulk GaAs (2023).
  * Banks  _et al._ [2017] H. B. Banks, Q. Wu, D. C. Valovcin, S. Mack, A. C. Gossard, L. Pfeiffer, R.-B. Liu, and M. S. Sherwin, Physical Review X 7, 041042 (2017).
  * Wu and Sherwin [2023] Q. Wu and M. S. Sherwin, Physical Review B 107, 174308 (2023).
  * Zaks  _et al._ [2013] B. Zaks, H. Banks, and M. S. Sherwin, Applied Physics Letters 102, 012104 (2013).
  * Borsch  _et al._ [2020] M. Borsch, C. P. Schmid, L. Weigl, S. Schlauderer, N. Hofmann, C. Lange, J. T. Steiner, S. W. Koch, R. Huber, and M. Kira, Science 370, 1204 (2020).
  * Langer  _et al._ [2016] F. Langer, M. Hohenleutner, C. P. Schmid, C. Poellmann, P. Nagler, T. Korn, C. Schüller, M. S. Sherwin, U. Huttner, J. T. Steiner, S. W. Koch, M. Kira, and R. Huber, Nature 533, 225 (2016).
  * Langer  _et al._ [2018] F. Langer, C. P. Schmid, S. Schlauderer, M. Gmitra, J. Fabian, P. Nagler, C. Schüller, T. Korn, P. G. Hawkins, J. T. Steiner, U. Huttner, S. W. Koch, M. Kira, and R. Huber, Nature 557, 76 (2018).
  * Banks  _et al._ [2013] H. Banks, B. Zaks, F. Yang, S. Mack, A. C. Gossard, R. Liu, and M. S. Sherwin, Physical Review Letters 111, 267402 (2013).
  * Hwang [1973] C. J. Hwang, Physical Review B 8, 646 (1973).
  * Becker  _et al._ [1988] P. C. Becker, H. L. Fragnito, C. H. B. Cruz, R. L. Fork, J. E. Cunningham, J. E. Henry, and C. V. Shank, Physical Review Letters 61, 1647 (1988).
  * Gopal  _et al._ [2000] A. V. Gopal, R. Kumar, A. S. Vengurlekar, A. Bosacchi, S. Franchi, and L. N. Pfeiffer, Journal of Applied Physics 87, 1858 (2000).
  * Rudin  _et al._ [1990] S. Rudin, T. L. Reinecke, and B. Segall, Physical Review B 42, 11218 (1990).
  * Alperovich  _et al._ [1976] V. L. Alperovich, V. M. Zaletin, A. F. Kravchenko, and A. S. Terekhov, physica status solidi (b) 77, 465 (1976).
  * Wang  _et al._ [1993] H. Wang, K. Ferrio, D. G. Steel, Y. Z. Hu, R. Binder, and S. W. Koch, Physical Review Letters 71, 1261 (1993).
  * Scholz [1995] R. Scholz, Journal of Applied Physics 77, 3219 (1995).
  * Ramian [1992] G. Ramian, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 318, 225 (1992).
  * Valovcin  _et al._ [2018] D. C. Valovcin, H. B. Banks, S. Mack, A. C. Gossard, K. West, L. Pfeiffer, and M. S. Sherwin, Optics Express 26, 29807 (2018).
  * Luttinger and Kohn [1955] J. M. Luttinger and W. Kohn, Physical Review 97, 869 (1955).
  * Vurgaftman  _et al._ [2001] I. Vurgaftman, J. R. Meyer, and L. R. Ram-Mohan, Journal of Applied Physics 89, 5815 (2001).