\section{Introduction}
\label{intro}
In their 2009 paper ``The Unreasonable Effectiveness of Data'' \cite{halevy2009unreasonable}, Alon Halevy, Peter Norvig, and Fernando Perreira urge ML researchers to ``follow the data and see where it leads.'' Since then, we have at an unprecedented scale amassed data for the purpose of machine learning. Yet, data quality---including the quality of human-collected data---has been left behind.

While there is a substantial body of work focusing on the reliability of human raters when performing evaluations, there are few \cite{diaz2018addressing, luo2020detecting, goyal2022your,prabhakaran2021releasing} if any studies investigating how the characteristics of rater pools impact ratings. That is, we know little about annotators’ individual characteristics (such as nationality, gender, education, race) and how they might influence the way they label data. This matters because, in seeking to build fair and responsible AI systems, we should anticipate potential biases that may emerge as a result of differences across user populations, and evaluation data should represent a variety of populations in order to better reflect viewpoints among real-world stakeholders and build diversity aware ground truth datasets\cite{aroyo2015truth}.

Responding to the \emph{Data-Centric AI} call to study impacts of data on AI systems \cite{Sambasivan_et_al}, we present findings from a semi-experimental exploration of rater diversity and its influence on safety annotations of chat bot conversations. We report results from a large-scale rater diversity study performed on a sample of 990 conversations generated by humans conversing with a generative AI-chatbot. We collected safety labels from a pool of 96 raters (recruited from two rating platforms covering a range of socio-economic subgroups) using 24 safety questions. Rather than requesting the typical number of ratings (a single rating or three-to-five ratings per conversation), we collected 40 ratings per safety question and repeated the experiment after 4 months. We analyzed:

\begin{itemize}[topsep=0pt,itemsep=0pt]
    \item the variance, disagreement and reliability within and across different demographic slices.
    \item the majority and minority perspectives on the data items.
    \item the consistency of raters across replications of the experiment.
\end{itemize}

The contribution and implications of our research are threefold. First, we propose a methodological approach for assessing the influences of rater diversity on data labels. Second, we demonstrate that raters' backgrounds do have an impact on the labels they provide. Third, we illustrate that there are discrepancies in rater consistency across two replications. As such, we contribute to two topics in this workshop: (1) repeatability of human evaluations and (2) methods for assessing the quality and the reliability of human evaluations.

\section{Experimental Design for Diverse Human Evaluation}
\label{experiment}
The overall research question guiding our study was: Do rater demographics help us improve reliability of annotations and ultimately provide guidance for collecting representative data for AI development and evaluation?

\paragraph{Data:}
\label{data}
The input data for this study was sampled from a 8K conversations dataset (composing 48K turns in total) generated by humans interacting with a generative AI-chatbot. In our study, the sample of 990 conversations was annotated for 24 safety characteristics (plus one general question about conversation understandability) by 40 unique raters per conversation, per safety characteristic. This annotation experiment was repeated with a different pool of 96 raters (including 19 overlapping raters in both replications) 4 months later. The total number of annotations gathered is close to 2 million (990K in each repetition).

\paragraph{Task:}
\label{task}
The safety evaluation task used a UI designed for the safety characteristics as defined in Appendix. Each rater was asked to rate the final chatbot response in each conversation as \emph{Unsafe}, \emph{Safe} or \emph{Unsure}. These answers were given for each of the 24 safety questions, which were organized in three groups (Q2-Q4 below), plus one question (Q1) checking the overall comprehensibility of the conversation.
\begin{itemize}[topsep=0pt,itemsep=0pt]
    \item Q1: Evaluate the whole conversation
    \item Q2: Evaluate the last bot utterance for “Harmful Content”
    \item Q3: Evaluate the last bot utterance for “Content with Unfair Bias”
    \item Q4: Evaluate the last bot utterance for “Misinformation and Political References”
\end{itemize}

\paragraph{Data Collection:}
\label{repetition}
We collected the ratings in two phases with an interval of four months. In both phases, we recruited 96 unique raters from two rater pools. All raters performed the task independently and used the same annotation template. All raters were asked to complete an optional demographic survey (e.g. gender, ethnicity, education level, age group, and native language). All questions in the demographic survey gave raters the option to select "Prefer not to answer". We also collected data about the average annotation time per conversation and the total time each rater spent annotating.

\paragraph{Raters:}
\label{raters}
In Phase 1, the breakdown of raters was: 71 in Pool 1 (42 India, 29 US), 25 in Pool 2 (12 India, 13 US). For Phase 2, the break down was: 72 in Pool 1 (40 India, 32 US), 24 in Pool 2 (12 India, 12 US). 19 of the raters participated in both phases (5 from Pool 1, 14 from Pool 2; the 5 raters from Pool 1 are all in US, 6 out of the 14 from Pool 2 are in India and 8 in US; 9 identify as female and 10 as male). In this paper, we report results from the Phase 2, however we compare the two phases to measure consistency for the raters who participated in both.

\section{Results}
\label{results}

We present four key high-level observations from this study that contribute to our understanding of reliability of human evaluations, and its relation to diversity among raters.

\paragraph{Unreliability of gold labels:} The left side of Fig. \ref{fig:combined} shows the difference between the number of raters saying \emph{Unsafe} vs. \emph{Safe} for each of the 990 conversations in our data. For around a quarter of the conversations, the number of \emph{Unsafe} and \emph{Safe} responses per conversation are quite similar, i.e., between 15-25 votes on either side. If only 3 to 5 annotators were to rate each item, as is common practice among researchers and practitioners building annotated datasets, this level of observed disagreement in these conversations may easily be lost. This suggests that majority-based (or even `unanimous') gold labels may be unreliable for a significant portion of the data, if the replication per item is low. This is a critical issue, since many evaluation tasks, even related to sensitive topics such as online safety, use such majority-based gold annotations to measure rater and model performance.
\begin{figure}[h!]
    \centering
    \includegraphics[width=0.43\textwidth]{figures/gold_vp.png}
    \includegraphics[width=0.53\textwidth]{figures/disagreementsperraterbins.pdf}
    \caption{\small \textbf{Left side:} Conversations are arranged horizontally, ranked by the difference between \emph{Unsafe} and \emph{Safe} votes. The y-axis shows this difference. The red square points to roughly a quarter of the conversations with nearly equal numbers of \emph{Unsafe} and \emph{Safe} votes. \textbf{Right side:} Histogram of the number of times each of the 19 raters who rated items in both phases disagreed with themselves.}
    \label{fig:combined}
\end{figure}

\paragraph{Inadequate intra-rater consistency:}
We now test whether the 19 raters who were present in both phases were consistent in their ratings. For the subset of items each of them annotated in both phases, we measured the number of times they disagreed with themselves (i.e., between phases) at least once --- considering all 25 questions separately. The right side of Fig. \ref{fig:combined} shows the histogram of disagreements for these 19 raters. Eleven of these raters disagreed with themselves at least ten, and as many as 3,199, times. This is another concerning finding that suggests there are extraneous factors that may significantly influence the consistency in raters' responses across different sittings at different points in time.

\paragraph{Disparate within-group coherence across subgroups:}
Despite the issues with consistency and reliability, we observed significant patterns in rater behaviour within and across the various subgroups we considered. For this analysis, we modeled each annotator's response to a conversation as a 72-dimensional \textit{response vector} that captures the one-hot encoding of the \{UNSAFE, SAFE, UNSURE\} answers for each of the 24 safety questions (Q2-Q4). This allows us to calculate the pair-wise distance between the response vectors of two raters as a metric for how strongly they disagreed with one another on any particular conversation prompt.

Fig.~\ref{fig:hamming} shows the average \emph{hamming distance} between all pairs of response vectors for each conversation, averaged across raters within different subgroups of raters. We observe that the average within-group rating distances vary
substantially across groups.

Lower hamming distance between a subgroup and \emph{All raters} means that the subgroup is consistent within itself and different than all raters. The results show disparities in agreement along three demographics - between US and Indian raters, between Pool 1 and Pool 2 and between female and make raters. In particular, US male raters in Pool 1 behaved more similarly among themselves than any of the other groups studied.

\begin{figure}[h!]
    \centering
    \includegraphics[width=.97\textwidth]{figures/hamming.png}
    \caption{The hamming distance metric on gender, locale and pool slices}
    \label{fig:hamming}
\end{figure}

\paragraph{Cross-group differences between subgroups:}
The above analysis provides only a partial picture, one that captures within-group distances but says nothing about whether the rating behaviors of a certain group of raters is more likely to be similar to others in the same subgroup than to those outside the subgroup. For instance, low within-group distance suggests that a particular subgroup has a coherent perspective on the task. If two different subgroups along a diversity axis (say, gender) exhibit such high within-group coherence, but also have low cross-group distance, it suggests that this particular diversity axis may not have any
substantial influence in the context of this task. However, if two groups that have low within-group distance, also have high cross-group distance, it suggests that the diversity axis is a
substantial
differentiator for the task.

Fig.~\ref{fig:cosine} shows the within-group and cross-group distance along the locale (IN vs. US), pool (Pool 1 vs. Pool 2), and gender axes. The results show that while US raters produced significantly more similar ratings with other US raters, compared to IN raters, on average. In the case of gender, female raters produced ratings that are very similar to each other, and significantly dissimilar to the ratings produced by male raters. Moreover, there was not much variance in the average distance across different female rater pairs, whereas male raters exhibited high variance across pairs in how much they disagreed with one another. While we observe some difference between the Pool 1 and Pool 2, those differences are not statistically significant.

\begin{figure}
    \centering
    \begin{subfigure}[t]{0.3\textwidth}
        \vskip 0pt
        \includegraphics[width=\textwidth]{figures/hamming_locale.png}

    \end{subfigure}\hspace{5mm}
    \begin{subfigure}[t]{0.3\textwidth}
        \vskip 0pt
        \includegraphics[width=\textwidth]{figures/hamming_platform.png}

    \end{subfigure}\hspace{5mm}
    \begin{subfigure}[t]{0.3\textwidth}
        \vskip 0pt
        \includegraphics[width=\textwidth]{figures/hamming_gender.png}

    \end{subfigure}
    \caption{Average pairwise hamming distance across different locale, platform, gender slices. The number of pairs within each group is specified in parenthesis.}
    \label{fig:cosine}
\end{figure}

\begin{comment}
The results from this study contribute to our understanding of reliability of human evaluations in relation to diversity among raters. We present four key observations:
\begin{itemize}[topsep=0pt,itemsep=0pt]
    \item 25\% of conversations rated by the diverse pool of raters had near-equal YES and NO ratings
    \begin{itemize}
        \item Gold annotations for safety (as determined by a simple majority-vote) would disagree with nearly half of the raters on these conversations
    \end{itemize}
    \item Raters are not always consistent with themselves when repeating the same task
    \begin{itemize}
        \item 20\% of raters disagreed with themselves on 10\% of the question-item pairs; Some disagree much more than others
    \end{itemize}
    \item Different demographic groups have different notions of safety:
    \begin{itemize}
        \item US raters overall rate differently than IN raters
        \item Female raters give very different results than all other raters
    \end{itemize}
        \item Different rater pools have higher internal agreement:
    \begin{itemize}
        \item US raters in pool 1 have the highest within-group agreement
        \item IN raters in pool 2 have the lowest within-group agreement
        \item Raters in pool 2 have high within-group disagreement both for IN and US
    \end{itemize}
\end{itemize}

\paragraph{Gold Annotations} Considering that many evaluation tasks use \emph{gold annotations} as a way to measure rater performance, our first observation is quite critical in this space. As we see in Fig. \ref{fig:gold} at least a quarter of the conversations in our dataset result in quite ambiguous annotations for safety (i.e., the number of UNSAFE and SAFE responses per conversation are quite similar --- in this case we see between 15-25 votes on each side). If only 3-5 annotators were to rate each item, such subtle levels of disagreement would be much less apparent.

\begin{figure}[h!]
    \centering
    \includegraphics[width=.97\textwidth]{figures/gold.png}
    \caption{Distribution of UNSAFE and SAFE votes on all conversations; the red square points to a roughly a quarter of the conversations with nearly equal votes}
    \label{fig:gold}
\end{figure}

To understand better the disagreement between raters, not only in these extreme cases, as shown above, we performed two similarity analysis across raters from different locales and raters from different platforms. First, we measured the \emph{in-group similarity} of the annotators belonging to the same subgroup with respect to the way they responded to the safety questions. Second, we measured the \emph{cross-group similarity} of a given subgroup with respect to how similar their answers were to outgroup raters.

\paragraph{Rater Consistency} Among the raters who annotated the same items in each replication, we measured the number of times they disagreed with themselves (i.e., between replications) at least once. Fig. \ref{fig:disbyrater} shows the histogram of disagreements for these 19 raters. Figure \ref{fig:disbyrater} shows the 19 annotators (about 20\% of all raters) who rated the same items in both replications. Eleven of these raters disagreed with themselves at least ten, and as many as 2,880, times.

\begin{wrapfigure}{r}{0.35\textwidth}
    \centering
    \includegraphics[width=0.35\textwidth]{figures/disagreementsperrater.pdf}
    \caption{Histogram showing the number of times each of the 19 raters who rated items in both replications disagreed with themselves.}
    \label{fig:disbyrater}
\end{wrapfigure}

\paragraph{Rater Similarity} To assess the similarity of two raters’ responses to the same conversation, we generated a 72-dimension \emph{response vector} that captures the one-hot encodings of the \{YES, NO, UNSURE\} answers for each of the 24 safety questions (i.e., Q2 onward) for each of the 990*40 annotations. Then, for each of the 990 conversations we calculated \emph{the cosine similarity} (see Fig. \ref{fig:cosine}) and \emph{hamming distance} (see Fig. \ref{fig:hamming}) between each pair of response vectors to that item, and average those scores across conversations. We also calculated the p-score measuring the likelihood of obtaining a more extreme score within a group under the assumption that there is no correlation between annotator response and demographics, using permutation tests. The results show significant disparities in agreement between US and Indian raters, and that both US male and female raters have higher agreement than the entire rater population. In particular, pool 1 annotators behaved more similarly than any of the other groups studied.

\begin{figure}[h!]
    \centering
    \includegraphics[width=.97\textwidth]{figures/hamming.png}
    \caption{The hamming distance metric on gender, locale and pool slices}
    \label{fig:hamming}
\end{figure}

\begin{figure}[b!]
    \centering
    \includegraphics[width=.97\textwidth]{figures/cosine.png}
    \caption{Cosine similarity distance on locale and pool slices}
    \label{fig:cosine}
\end{figure}

A cosine similarity distance of 0 between two items means they are different, and cosine similarity distance of 1 means they are identical. In \ref{fig:cosine}, all the cosine similarity distances for all slices are quite low (all less than 0.3), which indicates that raters from different locales and pools rate the same items quite differently. This supports the conclusion that even coarse demographic categories are quite meaningful to consider when evaluating results from human annotation.

\end{comment}

\section{Reflections}
\label{discussion}

In this paper we are excited to share just a few of the high-level results from the presented work. These results offer a clear indication that raters' demographics and the pool from which raters have been recruited have an impact on labelling tasks. Because the analysis we have done thus far is relatively coarse-grained, we believe that slicing further into the ethnicity, native languages and age groups of the raters is likely to reveal further insights and provide additional evidence of systematic differences between different groupings of raters. We will be conducting this detailed analysis with the ethnicity, age group and native language data that accompanies our data corpus and reporting results in upcoming publications.

As we propose a methodology for assessing the influences of rater diversity on data labels, our future work will also focus on determining the optimal number of raters per conversation and to what extent the impacts of rater diversity can be captured in smaller numbers of raters. This will be done in order to improve dataset generation methods that aim to address rater diversity.

Finally, we recognize more work is needed to help distinguish \emph{good} from \emph{bad} disagreement. In our work, this could be done by correlating the temporal data with other behavioral traits in raters across the two replications. Ultimately, this would extend our methodology to include an approach for studying outliers and different annotation perspectives.

\bibliographystyle{plain}


\section{Background}
\label{sec:intro}

During the last few decades, breast cancer has rapidly become a global healthcare burden, with an estimated 2.26 million new incidences and 0.68 million deaths in 2020 \cite{globocan}. Early detection and diagnosis of abnormalities in the breast might result in a better prognosis and eventually provide clinicians with better treatment options. In this regard, mammography is the gold standard breast cancer screening technology for women. Breast Imaging Reporting and Data System (BIRADS) is a tool used by radiologists to standardize and manage breast mammography screening examinations \cite{birads}. Abnormal findings such as mass and calcification are often regarded as high-risk indications and can reveal signs of malignancy, depending on their size, growth, and distribution. Based on the BIRADS score, a biopsy is recommended, where a small part of the tissue in the region of abnormality is removed and further analyzed using histopathological analysis. The pathologists then analyze the tissue using microscopy to confirm whether the finding is benign or malignant and if malignant, then further classify the finding into various pathological subtypes. Upon confirmation of malignancy, immunohistochemistry analysis is performed to investigate the molecular subtypes for finding the genetic expression of the cancer cells. The four main molecular subtype categories are luminal A, luminal B, human epidermal growth factor receptor 2 (HER2)-enriched, and triple-negative \cite {horvath2021molecular, cho2016molecular, boisserie2013correlation}. Finally, the clinicians use the molecular subtypes to prepare a personalized cancer treatment plan. Even though magnetic resonance imaging includes more information about the molecular subtypes, recent studies have shown imaging feature correlations between mammography imaging and molecular subtypes \cite {horvath2021molecular, cho2016molecular, boisserie2013correlation}.
Deep learning (DL)-based computer-aided diagnosis (CAD) solutions have been shown to achieve good performances in medical image analysis and, in many cases, have achieved on-par performance with clinicians \cite{nagendran2020artificial}. Most of the tasks involving DL methods in mammography analysis revolve around the classification, segmentation, and detection of abnormalities such as calcification, mass, and their benignancy and malignancy. However, there is a lack of research involving the analysis of molecular subtype features in full mammogram images using DL-based methods. A CAD system that can detect these molecular signatures directly from mammography can potentially streamline the clinicians' workflow, especially in countries where it is the only standalone cancer diagnosis tool. Furthermore, including additional information of these subtypes can result in designing a sophisticated multi-task learning-based CAD system.

\subsection{Previous Works}

Earlier works involving subtype analysis in full mammogram images mainly use handcrafted and radiomics features, where annotations from expert radiologists are acquired. These features are used to train machine learning models for classification \cite {son2020prediction, ma2022predicting, wu2019prediction}. Son, Jinwoo, et al. use a total of 129 radiomic features from 365 patients and classify the synthetic mammograms obtained from 3D breast tomosynthesis into luminal, HER2, and triple-negative subtypes \cite{son2020prediction}. Ma, Mengwei, et al. use the clinical information and 50 imaging features based on abnormal findings seen in mammograms and ultrasound imaging from 600 patients, to develop and compare different machine learning models for the subtype classification \cite{ma2022predicting}. The BIRADS information from MRI and mammography is converted into 82 imaging features from 363 patients and is used to develop a Decision Tree (DT) model by Wu, Mingxiang et al. \cite{wu2019prediction}. They show the combination of features extracted using both imaging modalities can boost the performance of the subtype classification model. By using BIRADS features, they indicate the generalization ability of the model, by being modality- and vendor-independent \cite{wu2019prediction}.

All the previous studies involving the subtype classification in mammography using deep learning methods are trained on the region of interest (ROI) from full mammogram images \cite {ueda2021training, muramatsu2022intrinsic, zhang2021predicting}. Ueda, Daiju, et al. proposed the first DL method to predict the hormone receptor expressions responsible for the molecular subtypes \cite{ueda2021training}. They use an ensemble of VGG, Inception Net, and ResNet models for training and validation of their method using 1448 ROIs from mammogram images \cite{ueda2021training}. Zhang, Tianyu et al. use an attention-mechanism-based ResNet-50 which uses a combination of mammogram and ultrasound ROIs as input to predict the four molecular subtypes. They also highlight the performance of luminal and non-luminal subtypes in their work \cite{zhang2021predicting}. Muramatsu, Chisako, et al. propose a contrastive pretraining on 385 ROI images, and k-nearest neighbor classifier on top of the extracted features for the classification of four molecular subtypes \cite{muramatsu2022intrinsic}. The datasets used in all the earlier works involve mammograms with high-quality pixel-level labels annotated by experienced radiologists and are not publicly available.

\subsection{Research Goals}
In this work, we investigate the feasibility of classifying the full mammogram images with only image-level labels into luminal and non-luminal subtypes using a convolutional neural network (CNN), trained and validated using the publicly available, Chinese Mammography Dataset (CMMD) \cite{cui2021chinese}. The main contributions of our study are: 1. Analysis of the complexity in training a CNN model, for luminal and non-luminal subtype classification; 2. Selection of a suitable CNN model for transfer learning (TL) by comparing three different CNNs for abnormality classification tasks; 3. Selecting the multi-label, multi-class (MLMC) ResNet-18 (R18)-based model to classify the calcification, mass, benign and malignant findings for transfer learning; 4. Developing a luminal subtype classifier using the pre-trained model from the previous step. To the best of our knowledge, we show the first baseline results on the CMMD dataset for the luminal versus non-luminal subtype classification.

\section{Methods}
\label{sec:methods}
\subsection{Data}

As mentioned previously, we use the publicly available CMMD dataset from the ``The Cancer Imaging Archive" (TCIA) \cite{cui2021chinese, clark2013cancer}. The dataset consists of 3,744 mammograms belonging to 1,775 patients from China. Each mammogram is labeled with the age of the patient, craniocaudal (CC)/mediolateral oblique (MLO) views, right/left breast, calcification/mass anomalies, and benign/malignant lesions. A subset of these mammograms from 749 patients with malignant lesions is labeled with luminal A, luminal B, HER2-enriched, and triple-negative subtypes, respectively. In this work, we aim to perform a luminal versus non-luminal classification study. Therefore, luminal A and luminal B subtypes are combined to form the `luminal' class. HER2-enriched and triple-negative are merged into the `non-luminal' class. We use the entire dataset for the abnormality classification tasks and the subset for the luminal versus non-luminal classification task.

\begin{table}[!ht]
    \centering
    \caption{Class sample distribution in the CMMD dataset. The entire dataset is split into approximately 90:10 for training and testing. 90\% of the dataset is used for the 5-fold cross-validation (CV) splits, where almost 80\% is used for training and 20\% is used for validation.}
    \label{tab:dataset}
    \begin{subtable}{\linewidth}
    \centering
        \begin{tabular}{c  c  c }
         \toprule
                       & 5-fold CV Dataset          &  Test Dataset \\
        \midrule
         Luminal         & 956                        & 100            \\
        Non-Luminal         & 404                        & 38            \\
        \midrule
       Total           &  1360                      & 138         \\

         \bottomrule
        \end{tabular}
        \caption{Dataset used for luminal versus non-luminal classification task.}
        \label{tab:dataseta}
    \end{subtable}

    \begin{subtable}{\linewidth}
        \centering
        \begin{tabular}{c  c  c }

        \toprule
                & 5-fold CV Dataset        &  Test Dataset \\  [0.5ex]
         \midrule
         Benign          &  1004                    & 108           \\

        Malignant       &  2386                    & 246           \\

        \midrule

        Calcification   &  476                     &  48                    \\

        Mass            &  2076                    & 222                 \\

        Calc+Mass       &  838                     & 84                 \\

        \midrule
        Total           &  3390                    & 354                \\

         \bottomrule

        \end{tabular}
        \caption{Dataset used for the abnormality classification tasks.}

        \label{tab:datasetb}
    \end{subtable}
\end{table}

\cref{tab:dataset} shows the distribution of the class samples in both the datasets used in this study. We perform a 5-fold cross-validation (CV) in all experiments conducted in this study and test individual folds using the independent test dataset. 90 \% of the total dataset is used for CV and 10 \% is used as the independent test dataset. Approximately, 80\% of the CV dataset is used for training, and 20 \% is used for validation across all the folds. We ensure that all the folds consist of an almost equal distribution of the class samples and that there is no leakage of data samples belonging to the same patient between the training and validation sets. The mammograms in the CMMD dataset are converted from Digital Imaging and Communications in Medicine (DICOM) files into JPEG format. The quality of the saved images is preserved by scaling the pixel values using the window-width and window-level from the DICOM metadata. Furthermore, the black pixels outside the breast region are cropped-out to obtain only the breast view as the input to the CNN. The images are resized to 1326 x 512 using the average aspect ratio of all cropped images in the dataset to meet the computational requirements for model training.

\begin{figure*}[htb]

\begin{minipage}[b]{1.0\linewidth}
  \centering
  \centerline{\includegraphics[width=11.7cm, height=6.7cm]{isbi_mlmc.pdf}}

\end{minipage}
\caption{Overview of the transfer learning method: MLMC classification network is used for supervised pretraining for the mass, calcification, benign, and malignant classes. The pretrained network is then used for finetuning the CNN for the luminal subtype classification task.}
\label{fig1}
\end{figure*}

\subsection{Experimental Setup}

R18 pretrained on ImageNet is used as the standard CNN architecture across all the experiments conducted in this research. The CNN is selected after searching for the right balance between the number of training parameters and FLOPs to ensure faster training performance on large-resolution input images. During training, we perform four data augmentations, namely random horizontal flip (p=0.5), AugMix (p=0.2) \cite{hendrycks2019augmix}, random histogram equalization (p=0.4), and random erasing (p=0.1). We employ an early stopping criterion with a patience value of 10, to prevent the models from overfitting and achieve faster training times. All models are saved on the lowest validation loss to avoid using the overfitted models during testing. The experimental setup in this study is designed using the PyTorch framework.

\subsubsection{Baseline: Luminal Versus Non-Luminal Classification Task}

\label {sub:baseline}

At first, we perform the luminal versus non-luminal subtype classification, by refactoring the fully connected (FC) layer of R18 with two output nodes, representing the two classes. As the training is extremely unstable, we perform experiments by adding the dropout layer with various probabilities, eventually selecting a value of 0.3. A weighted random sampler is used during training to address the data imbalance problem in the subtype dataset. Based on our initial experiments, we see a significant gain in performance while using this sampling strategy. A batch size of 16 is used for all the experiments and a cross-entropy (CE) loss is optimized using the standard Adam optimization strategy with a learning rate of $1e^ {-5} $ and a weight decay of $5e^ {-3} $. It is seen in our experiments that using low learning rates can stabilize the training to some extent for this task.

\subsubsection {Abnormality Classification Tasks: Binary-Class And Multi-Label Multi-Class Task}

As the performance of CNN is poor on the test dataset and training instabilities are observed during the 5-fold CV on the first task, we investigate whether TL from a related task can boost the performance of the luminal classifier. For this purpose, we investigate three classifiers, namely, mass versus calcification, benign versus malignant, and a combined MLMC classifier that classifies mass, calcification, and its benignancy and malignancy. The goal here is to perform supervised training on the aforementioned tasks and select the right model to perform TL. As the imaging features correlated to the subtypes arise due to abnormalities such as mass and calcification \cite {horvath2021molecular, cho2016molecular, boisserie2013correlation}, we initially employ a mass versus calcification classifier. As many mammograms consist of both mass and calcifications, a multi-label binary classifier is used with a binary CE loss function. The benign versus malignant classification is then explored, as the molecular subtypes are the hormone receptor expressions arising from the malignant lesions. We employ a binary output on top of the FC layer and use a CE loss function. Finally, the MLMC classifier is employed with the intuition that a single network should be able to learn robust representations of both the mass versus calcification and benign versus malignant task. In this case, the FC layer is refactored to output three neurons, each representing the probability of one of the three classes. As only a subset of the mammograms with malignant lesions is labeled with the molecular subtypes, we do not include the subtype label in the MLMC task. A binary CE loss is employed again in this case. Based on extensive experimentation as reported in \cref{sec:results}, we propose to use the MLMC classifier for supervised pretraining. A learning rate of $1e^ {-4} $ and weight decay of $5e^ {-3} $ is used for the Adam optimizer.

\subsubsection {Transfer Learning Task: Luminal Versus Non-Luminal Classification Task}

The overview of the approach of the TL task can be seen in \cref{fig1}. We replace the fully connected layer from three to two, to match the luminal and non-luminal subtype classes. Each fold is initialized by the pre-trained weights from each fold in the MLMC task, and all layers in the R18 model are finetuned to ensure optimal performance. The remaining settings from \cref{sub:baseline} are retained for a fair comparison between the two methods.

\section{Results and Discussion}
\label{sec:results}

\begin{figure*}[htb]

\begin{minipage}[b]{1.0\linewidth}
  \centering
  \centerline{\includegraphics[width=13.2cm, height=6cm]{ISBI_GradCAM.pdf}}
\end{minipage}
\caption{Grad-CAM visualizations and magnified ROI of the CC and MLO view from the mammogram of a 58-year-old breast cancer patient with both calcification and mass findings. The malignant lesion belongs to the luminal subtype.}
\label{fig:grad}
\end{figure*}
The per-class F1 score, total F1 score, and Area Under the Curve (AUC) score are used as the evaluation metrics for the luminal subtype classification. As there are three separate tasks in abnormality classification, only the per-class F1 score is used to evaluate the performance of the models. The model from each of the five folds is used for testing and the mean and standard deviation values are reported.

From \cref{tab:restable2}, it is clear that the luminal versus non-luminal classification model fails to learn good representations of the molecular subtypes from the full mammogram images. It achieves a mean F1 score of 0.4987 and a mean AUC of 0.5358, with a mean F1 score as low as 0.2099 for the non-luminal class. One potential reason for the inferior performance can be attributed to the class imbalance, especially with the high F1 score of the dominant luminal subtype class, even after using weighted random sampling. Moreover, there is no additional clinical information on the local appearance of calcification and mass in the CMMD dataset. For example, the presence of microcalcification clusters, mass with well-circumscribed or poorly-circumscribed margins, and spiculations are some of the imaging features which highlight the specific molecular subtype signatures \cite{cho2016molecular}. As this combination of clinical knowledge in the form of class label information is missing, the network fails to learn the correlations between the abnormalities and molecular subtypes. This can also be supported by the increased performance of the CNNs as reported in the literature when trained on smaller ROIs~\cite{ueda2021training, zhang2021predicting}.

\begin{table}[t]
    \centering
    \caption{Classification performance on the test dataset. The mean values on the test dataset from 5-fold CV are reported. Standard deviations are reported in brackets. The best results are shown in bold.}
    \begin{subtable}{\linewidth}
        \centering
        \begin{tabular}{c  c  c }
         \toprule
        Class                & Binary-class F1         & MLMC F1 \\
        \midrule
         Calcification       & 0.7765 (0.0116)                   & \textbf{0.8025     (0.0190)}           \\

        Mass                 & 0.9233  (0.0036)                   & \textbf{0.9238   (0.0039)}    \\

        Benign               & 0.4079  (0.083)                  &  \textbf{0.5061    (0.0860)}                 \\

        Malignant           & 0.8313  (0.0039)                   &  \textbf{0.8480   (0.0079)}          \\
         \bottomrule
        \end{tabular}
        \caption{Binary versus MLMC results.}
        \label{tab:binmlmc}
    \end{subtable}\\%
    \begin{subtable}{\linewidth}
        \centering
        \begin{tabular}{c  c  c }

        \toprule
                                & Baseline               &  Transfer learning \\ [0.5ex]
         \midrule
        Luminal F1              & 0.7864 (0.0667)      & \textbf{0.8188     (0.0245)}       \\

        Non Luminal F1          & 0.2099 (0.1749)     & \textbf{0.5199 (0.0253)}       \\

        F1 Score                & 0.4987 (0.0647)      &  \textbf{0.6693  (0.0184)}       \\

        AUC                     & 0.5358 (0.0300)      & \textbf{0.6688     (0.0210)}    \\
         \bottomrule

        \end{tabular}
        \caption{Baseline versus Transfer learning results.}
        \label{tab:restable2}
    \end{subtable}
\end{table}

In comparison, it is easier for the network to classify the calcification, mass, benign, and malignant features. This can be seen quantitatively in \cref{tab:binmlmc} and might be due to the global appearance of calcification and mass on mammograms. Calcifications appear as bright, high-intensity, white spots on mammograms, while the global morphological features of most of the mass abnormalities are distinct. Therefore, we can hypothesize that these abnormalities are representative of their respective class labels. The same applies to malignant and benign classes as well. Moreover, we observe the network works better when we train a single network to classify all the abnormalities. Although the performance of the mass classification stays the same, we obtain more than 2\%, 9\%, and 1\% performance boost for the calcification, benign and malignant classes, respectively. We find these gains in performance to be statistically significant while using a t-test, with a two-tailed $p$-value of $p$$<$0.0001. Hence, we use the MLMC model for finetuning the luminal classifier.

We achieve a mean AUC of 0.6688 and a mean F1-score of 0.6693 in the TL task. A high-performance boost with statistical significance ($p$$<$0.0001) in comparison to the baseline is obtained by performing TL, and this trend is seen across all the 5-folds. This can be observed in \cref{tab:restable2}, where we see that the increase in the overall performance can be attributed to the gains in the non-luminal class, showing the network no longer favors the majority class. The clinical information in the form of missing additional labels in the baseline models is replaced by rich representations of the same information learned from the supervised pretraining in the MLMC task.
We further analyze the Gradient-Weighted Class Activation Mappings (Grad-CAM) to visualize the attention of the MLMC network. \cref{fig:grad} shows an example of the Grad-CAM obtained from the MLMC model. The suspected ROI is comparably found by the Grad-CAM in both CC and MLO views. Even though the network predictions are the same as the input class label, and the ROI looks visually suspicious by clinical definitions of mass and calcification, it is difficult to use these as weak localization maps for further analysis due to the absence of lesion location data.

\section{Conclusion and Future Work}
\label{sec:typestyle}
In this work, we investigate the luminal versus non-luminal subtype classification on full mammogram images using only image-level labels. We provide initial insights into the issues with CNNs for this task and also show that TL from a robust MLMC abnormality classification model can significantly boost the performance of the subtype classification. We achieve a mean AUC and F1 score of 0.6688 and 0.6693 using TL in comparison to 0.5358 and 0.4987 while not using TL, respectively. These initial results on the luminal classification tasks are promising and reiterate the need for more active research in this sub-area of breast cancer to eventually develop a standalone multi-task CNN model for analysis of all types of breast abnormalities. In the future, we plan to validate the Grad-CAM localization results from a subset of this data with the help of a clinician and use it as the test dataset. We then aim to use the Grad-CAM generated localization maps to develop a weakly-supervised ROI classifier model and also extract radiomics features for further comparison with the current methods in luminal subtype analysis.

\section{Compliance with ethical standards}
\label{sec:ethics}

This research study is conducted using the open access data available in ``The Cancer Imaging Archive" (TCIA) \cite{cui2021chinese, clark2013cancer} and consists of human subject data.
The license associated with the open-access data confirms the non-requirement of ethical approval.

\bibliographystyle{IEEEbib}


\section{Introduction}

Optical frequency combs (OFCs) denote light fields consisting of a discrete series of evenly spaced spectral components (or colors) that maintain high coherence across the whole bandwidth with ultra-low-noise~\cite{CY-RMP-03}. Their applications encompass a variety of fields including ultra-broadband coherent optical communications, high-precision optical spectroscopy, high precision metrology and optical arbitrary pulse generation~\cite{UHH-NAT-02,CW-NAP-10,D-JOSAB-10,PPR-PR-18}, to name just a few. Prominent examples of OFCs realizations leverage the temporal cavity solitons obtained in injected passive high-Q Kerr resonators, such as microrings \cite{HBJ-NAP-14} and fiber loops \cite{LCK-NAP-10} or employ active optical sources such as mode-locked lasers~\cite{tropper04,LMB-OE-10,WVM-LSP-17}.

Recently, the promising potential of injected Kerr micro-cavities enclosed into a long external feedback cavity as sources of high-power tunable OFCs, temporal localized structures (TLSs) and square-waves (SWs) was demonstrated~\cite{SPV-OL-19,SJG-OL-22,KSG-OL-22}. There, it was shown that the natural modeling approach consists in using singularly perturbed time-delayed systems, in which the algebraic constraints simply corresponds to the boundary conditions linking the optical field in the various parts of the compound cavity system. In the normal dispersion regime and for the resonant optical feedback, dark and bright TLSs (and corresponding OFCs) appear via the locking of domain walls that connect the high and low bistable continuous wave (CW) intensity levels. Due to the oscillatory tails induced by the cavity dispersion~\cite{SCM-PRL-19}, these solutions can interlock at multiple equilibrium distances leading to a rich solution multistability. Moreover, operated out of the bistability regime and for anti-resonant feedback, SWs possessing a period of twice the round-trip in the external cavity were disclosed. Interestingly, the latter appear through a homoclinic snaking scenario leading to the formation of complex-shaped multistable solutions.

Departing from former studies, in this paper we analyze the impact of considering a semiconductor quantum well as the nonlinear element. A first principle model for the optical response is employed which allows to explore e.g., the influence of the detuning with respect to the semiconductor band-gap that controls the ratio between nonlinear index change and absorption. Using a combination of analytical, numerical and path-continuation methods, we shall show that this extended model predicts the existence of a multistable set of bright and dark TLSs coexisting on their respective bistable homogeneous CW backgrounds. Furthermore, in the case of anti-resonant optical feedback, the dynamics of resulting SWs is analyzed. Finally, in order to clarify the influence of the second and third order chromatic dispersion and of the frequency dependence of the quantum-well response, we perform a multiple time-scale analysis in the good cavity limit. The resulting dispersive normal form partial differential equation shows a good qualitative agreement with the original time-delayed model.

\section{Model equations}

\begin{figure}
	\includegraphics[width=1\columnwidth]{fig1.pdf}
	\caption{a) Schematic of a micro-cavity containing a quantum well layer coupled to an external cavity with round trip $\tau$ which is closed by a
mirror with reflectivity $\eta$ and phase $\phi$, and driven by a CW
beam with amplitude $Y_0$. b) A single TLS circulating in the external cavity obtained by a numerical integration of Eqs.~\eqref{eq:QGTI1}-\eqref{eq:QGTI3}. The inset shows the dynamics over several round-trips. c) Power density spectrum of the solution in b). Parameters are $(\delta,h,\eta,\varphi,\gamma,f,u,b,\tau)=(0.015,2,0.99,0,1,0.2,0.65,1000,250)$.} \label{fig:fig1}
	\end{figure}

The schematic setup is depicted in Fig.~\ref{fig:fig1}~(a). It is composed of a monomode micro-cavity of a few micrometers in length, with a radius up to $100\,\mathrm{\mu m}$ and with round-trip time $\tau_c$, which contains a thin quantum well layer acting as a nonlinear medium that is situated at the anti-node of the field. The micro-cavity is closed by two distributed Bragg mirrors with reflectivities $r_{1,2}$ and it is coupled
to a long external cavity of (typically) a few centimeters with
round-trip time $\tau\gg\tau_c$. The system is subjected to CW injection with
amplitude $Y_0$ and frequency $\omega_0$. The detuning with respect to the
micro-cavity resonance $\omega_c$ is $\delta=\omega_c-\omega_0$. The external cavity is closed by a feedback mirror with reflectivity $\eta$. The total external cavity phase $\varphi$ consists of the accumulated phase per round-trip due to propagation and the phase of the feedback mirror $\phi$.

Following the methods developed in Refs.~\cite{MB-JQE-05,SPV-OL-19}, the first-principles model for the system shown in Fig.~\ref{fig:fig1}~(a) reads
\begin{align}
        \dot{E}&=\left[if\chi -1 -i\delta\right]E+hY\,, \label{eq:QGTI1}\\
        \dot{D}&=\gamma\left[-D+\Im(\chi)|E|^2\right]\,, \label{eq:QGTI2}\\
        Y&=\eta e^{i\varphi}\left[E(t-\tau)-Y(t-\tau)\right]+Y_0\sqrt{1-\eta^2}\,.\label{eq:QGTI3}
\end{align}
Here, $E$ and $Y$ denote the slowly varying envelope of electric
fields in the micro-cavity and the external cavity, respectively. Further, $D$ is the carrier density normalized to the transparency level at the band gap, $\gamma$ is the scaled carrier lifetime, and $f$ is an effective dimensionless gain prefactor. The latter is the product of a geometrical factor that contains the cavity enhancement effect and the active medium length, normalized by the photon lifetime in the micro-cavity.
The field cavity enhancement can be conveniently scaled out using the Stokes relations allowing the fields $E$ and $Y$ to be of the same order of magnitude. This leads to a simple input-output relation $O = E - Y$: The output $O$ is the combination of the intra-cavity photons transmitted and reflected by the micro-cavity and it is re-injected after attenuation and a time delay $\tau$. The coupling between the intra- and external cavity fields is given by Eq.~\eqref{eq:QGTI3} which is a delay algebraic equation (DAE). The DAE~\eqref{eq:QGTI3} takes into account all the multiple reflections in a possibly high finesse external cavity where $\eta\lesssim 1$.

The light coupling efficiency in the cavity is given by the factor $h=h(r_1,\,r_2)=(1+|r_2|)(1-|r_1|)/(1-|r_1||r_2|)$ in Eq.~\eqref{eq:QGTI1}. Here, we consider the case of a perfectly reflecting bottom mirror, i.e. $h=h(r_1,\,1)=2$ which corresponds to the so-called imbalanced Gires–Tournois interferometer regime~\cite{GT-CRA-64}. These interferometers are widely used as optical pulse shaping elements~\cite{GAS-OL-00}. They induce a controllable amount and sign of group delay dispersion (GDD), which is typically the dominating effect outside resonance. The amount and sign of GDD are tunable by choosing the frequency of operation with respect to the cavity resonance and around resonance where the second-order contribution vanishes and switches sign, third order dispersion (TOD) becomes the leading term~\cite{SCM-PRL-19}. Due to TOD, the resulting temporal pulses can possess oscillatory tails, see Fig.~\ref{fig:fig1}~(b), where a typical periodic time trace obtained from integrating numerically Eqs.~\eqref{eq:QGTI1}-\eqref{eq:QGTI3}  is shown. As a result, the envelope of the corresponding OFC is asymmetrical, see Fig.~\ref{fig:fig1}~(c).

Finally, $\chi=\chi(u,D)$ is a nonlinear susceptibility function. As a simple model, we employ the quantum well susceptibility as introduced in Ref.~\cite{B-PRA-98} as
\begin{align}
        \chi=\frac{1}{\pi}\left[\ln(u+i-b)-2\ln(u+i-D)+\ln(u+i)-\ln(b)\right].\label{eq:QGTI_chi}
\end{align}
The susceptibility $\chi$ describes the response of the medium to incoming light with frequency $\omega_0$. This analytical formula originates from assuming two parabolic bands with a quasi-equilibrium Fermi distribution for electron and holes.
This response possesses both a real and an imaginary part that represents the non-instantaneous Kerr effect mediated by the population inversion, as well as the field saturable absorption, which also depends on $D$.
\begin{figure}
	\includegraphics[width=1\columnwidth]{fig2.pdf}
	\caption{The real (red) and negative imaginary (blue) parts of the
	susceptibility $\chi\left(u,D\right)$ as given by Eq.~\eqref{eq:QGTI_chi} for $b=100$, for $D=0$ (solid line) and $D=0.3$ (dash-dotted line). \label{fig:fig2}}
	\end{figure}
Further, $\chi$ depends on $u$, the normalized frequency with respect to the band-gap frequency $\omega_G$ as $u=(\omega_0-\omega_G)/\gamma_{\perp}$, where $\gamma_{\perp}$ is the intraband polarization dephasing rate. The parameter $b$ relates to the top-band frequency $\omega_T$, and reads $b=(\omega_T-\omega_G)/\gamma_{\perp}$. The gain (or the saturable absorption) is defined as $-\Im(\chi)$ while the phase shift caused by the nonlinear refraction index is described by $\Re(\chi)$.

The situation where $D=0$ describes an empty conduction band and a transition from transparency to saturable absorption located around $u=0$ while $D=1$ corresponds to the optical transparency being reached at the bandgap, i.e. $\Im[\chi(0,1)]=0$. The real and negative imaginary parts of the susceptibility $\chi$ are represented for different values of the population inversion in Fig.~\ref{fig:fig2}.

\section{Bifurcation analysis of the DAE system}\label{sec:DAE}

To understand the origin and the behaviour of the periodic pulse solution presented in Fig.~\ref{fig:fig1}~(b), we perform a bifurcation analysis of the system~\eqref{eq:QGTI1}-\eqref{eq:QGTI3} using a recently developed extension of DDE-BIFTOOL~\cite{DDEBT} that allows for the bifurcation analysis of algebraic and
neutral delayed equations~\cite{HGJ-OL-21,SJG-OL-22,KSG-OL-22,SGJ-PRL-22,BKW_JDEA_06}. We start the analysis with the situation where the injection is set to be resonant with an external cavity mode, i.e., $\varphi=0$. In the long-delay limit, where the time delay is the longest relevant time scale and in conjunction with sufficiently large detuning $\delta$, the quantum well nonlinearity causes a bistable CW response of the system in a certain parameter range of the injection intensity. We show in Fig.~\ref{fig:fig3}~(a) the resulting hysteresis diagram depicted in blue for the same parameter set used
in Figs.~\ref{fig:fig1}~(b,c).
The stability is indicated by the line thickness. One can see that the CW branch loses stability in a saddle-node (SN) bifurcation at the point $SN_1$ and regains it in another SN bifurcation at $SN_2$. From the CW branch, one can see that the branch of the TLSs emerges in a subcritical Adronov–Hopf bifurcation. Indeed, in the bistable CW region, similar to the case of the injected microcavity subjected to the Kerr nonlinearity~\cite{SPV-OL-19,SJG-OL-22}, the coexistence of two stable CW solutions can lead to the formation of domain walls between them that exhibit oscillatory tails, for certain sets of parameters. In particular, two opposed fronts may lock at discrete distances forming stable bound states which results in the formation of TLSs. The corresponding branch of periodic solutions (cf. red line in Fig.~\ref{fig:fig3}~(a)) features stable dark and bright TLSs over an extended injection range. Here, different stable and unstable locking ranges occur along a collapsed snaking structure~\cite{KW_PD_05,PKG_PRA_16} around a Maxwell line where opposing fronts have the same drifting speed, see the insets in panel (a). On the left or the right of the Maxwell line, the drifting speeds differ, leading to a transient towards the lower or the upper CW state, respectively. Figure~\ref{fig:fig3}~(b,c) show two exemplary dark and bright solution profiles at different positions along the snaking branch. Remarkably, one notices that although being calculated for very close values of $Y_0$, the dark TLS shows a stronger oscillatory tail, while the snaking is more pronounced in the vicinity of the lower CW state, i.e., for the bright solutions.
  \begin{figure}
	\includegraphics[width=1\columnwidth]{fig3.pdf}
		        \caption{(a) Branches of CW (blue) and TLSs (red) solutions of Eqs.~\eqref{eq:QGTI1}-\eqref{eq:QGTI3} as a function of the injection $Y_0$. Thick (thin) lines denote stable (unstable) solutions, respectively. The insets show zooms around the upper and lower snaking regions. In the right panels two exemplary dark (b) and bright (c) solution profiles at positions labeled with a red and green star in (a) are displayed. Parameters as in Fig.~\ref{fig:fig1}~(b,c).\label{fig:fig3}}
	  \end{figure}

To investigate the existence and the stability borders of the obtained TLSs solutions in this parameter regime, four two-parameter continuations in $(Y_0,\,\delta)$, $(Y_0,\,f)$, $(Y_0,\,u)$ and $(Y_0,\,\eta)$ planes are carried out and the results are presented in Fig.~\ref{fig:fig4} (a)-(d), respectively.
\begin{figure}
	\includegraphics[width=1\columnwidth]{fig4.pdf}
		        \caption{(a)-(d) Two-parameter bifurcation diagrams in  $(Y_0,\,\delta)$, $(Y_0,\,f)$, $(Y_0,\,u)$ and $(Y_0,\,\eta)$ planes, respectively. Evolution of the points $SN_1$ and $SN_2$ limiting the region of CW bistability (black) is shown together with the two major folds limiting the stability of the dark (dark red) and bright (light blue) TLSs, respectively. The horizontal dashed lines indicate the particular parameter values used in Fig.~\ref{fig:fig3}. Other parameters as in Figs.~\ref{fig:fig1},~\ref{fig:fig3}.\label{fig:fig4}}
	  \end{figure}
Here, black lines correspond to the thresholds of the SN bifurcation points $SN_1$ and $SN_2$ limiting the region of CW bistability, whereas dark red and light blue lines show the evolution of the two major folds limiting the stability of the dark and bright TLSs, respectively, cf. $F_{1,2}^{t,b}$ in the insets of Fig.~\ref{fig:fig3}~(a). The horizontal dashed line indicate the particular parameter values used in Fig.~\ref{fig:fig3}.

In particular, in Fig.~\ref{fig:fig4}~(a) one can see that the CW bistabile region exists for both anomalous and normal dispersion regimes and the stability regions of both bright and dark TLSs is broader for small detunings and injection values. The latter vanish for large values of $\delta$ and $Y_0$. A similar trend is presented in panel (c), where the diagram in the $(Y_0,\, u)$ plane is shown. Here, one can see that the bistability region opens only in a certain range around $u\sim 1$, which corresponds to the most efficient nonlinear modulation of the losses. Again, the region of stable existence for TLSs decreases for increasing injection values. Figure~\ref{fig:fig4}~(b) presents the evolution of the corresponding folds in the effective gain prefactor $f$. Here, one can see that the bistability region starts to exist for small finite values of $f$ and small detuning values and it is becoming broader with increasing $f$, hence for stronger nonlinearty; the fold lines of both dark and bright TLSs follow the same trend. In panel (d), the dynamics of the limiting folds is shown in the $(Y_0,\,\eta)$ plane. One can see that the TLSs  solutions exist and are stable only in the so-called good cavity limit, where $\eta\rightarrow 1$. Moreover, the corresponding regions are moving in the direction of larger injections with increasing $\eta$. Note that this analysis was carried out for $\gamma=1$, i.e. assuming that the carrier lifetime is equal to the photon lifetime in the cavity, a situation that typically corresponds to a fast saturable absorber with a relaxation rate in the picosecond range.

Finally, we present in Fig.~\ref{fig:fig5} a continuation in the carrier lifetime $\gamma$. Here, the evolution of the major folds limiting the stability of the dark (dark red) and bright (light blue) TLSs is shown, cf. the inset in Fig.~\ref{fig:fig3}~(a). One can see that for increasing $\gamma$ values, first bright TLSs start to exist and the bistability between dark and bright TLSs is established only above some finite value of $\gamma$. Moreover, one can see that the width of the stability regions remains constant for large $\gamma$ values, where the adiabatic elimination of the carriers can be performed, connecting our results with past studies based upon an instantaneous nonlinearity~\cite{SPV-OL-19,SJG-OL-22,KSG-OL-22}.
 \begin{figure}
	\includegraphics[width=1\columnwidth]{fig5.pdf}
		        \caption{Two-parameter bifurcation diagram in the $(Y_0,\,\gamma)$ plane showing the evolution of the major folds of the dark (dark red) and bright (light blue) TLSs. The horizontal black line indicates the particular parameter value used in Fig.~\ref{fig:fig3}. Other parameters as in Figs.~\ref{fig:fig1},~\ref{fig:fig3}\label{fig:fig5}.}
	  \end{figure}

\section{Dispersive master equation}\label{sec:PDE}
Although the bifurcation analysis allows for quantitative predictions of the behavior of the DAE system~\eqref{eq:QGTI1}-\eqref{eq:QGTI3}, an intuitive interpretation of the dynamics presented in Fig.~\ref{fig:fig3} is difficult since the role and impact of e.g., GDD and TOD on the TLS formation is hidden. Hence, to obtain a better understanding of the underlying physical mechanisms, we derive a partial differential equation (PDE) that approximates the dynamics of the full DAE model~\eqref{eq:QGTI1}-\eqref{eq:QGTI3} in the limit of a good cavity ($\eta\rightarrow 1$) and small detunings $\delta$ using a rigorous multiple time scale analysis as discussed in~\cite{SCM-PRL-19,SHJ-PRAp-20,SGJ-PRL-22,SJG-OL-22}.
We start our analysis by eliminating the field $Y(t)$  and transforming the DAE to a neutral delay differential equation (NDDE) that reads
\begin{eqnarray}
\eta\dot{E}(t-\tau)&+&\dot{E}(t)=\eta\left[if\chi(D(t-\tau))-1-i\delta+h\right]E(t-\tau)\nonumber\\
&+&\left[if\chi(D(t))-1-i\delta\right]E(t)+hY_0\sqrt{1-\eta^2}\,. \label{eq:NDDE}
\end{eqnarray}
Further, rescaling time with respect to the time-delay $\tau$ reveals a natural smallness parameter $\epsilon=\tau^{-1}$, where $\epsilon\ll 1$ in the long-delay limit. In the next step rescalings are chosen as  $E\rightarrow\epsilon^\frac{1}{2}\mathcal{E}$, $Y\rightarrow\epsilon^\frac{1}{2}\mathcal{Y}$, $Y_0\rightarrow\epsilon^\frac{1}{2}\mathcal{Y}_0$, and $D\rightarrow\epsilon\mathcal{D}$.  This choice of the scaling is also reflected in the treatment of the nonlinearity which now becomes $\chi(\epsilon\mathcal{D})$. Hence, the obvious next step is an expansion of $\chi$, leading to
\begin{equation}
 \chi(\epsilon\mathcal{D})=\chi(0)+\epsilon\chi^\prime(0)\mathcal{D}+\mathcal{O}(\epsilon^2)\,,
\end{equation}
where we ommited the $u$ dependence for brievity and
\begin{align}
\chi(0)=\chi_0&=-\frac{1}{\pi}\ln(-u-i)\nonumber\\
&=\frac{-1}{\pi}\ln\sqrt{u^2+1}+i\left[\frac{1}{2}+\frac{1}{\pi}\arctan(u)\right]\,,\\
\chi^\prime(0)=\chi_0^\prime&=\frac{2}{\pi}\frac{u-i}{u^2+1}\,.\nonumber
\end{align}
The dynamics is now considered on multiple timescales in powers of $\epsilon$.
 To this aim, the remaining fields and parameters are all expanded in powers of $\epsilon$. In addition, the fastest timescale $\sigma=\omega t$ is also stretched by expanding $\omega$ in powers of $\epsilon$ in order to account for the fact that the period $T$ is slightly larger than the time delay $\tau$. Hence, $\sigma$ describes the fast dynamics within one period.
 %
 \begin{figure}
\includegraphics[width=1\columnwidth]{fig6.pdf}
\caption{\label{fig:test1}  CW solutions of the PDE system~\eqref{eq:PDE1}-\eqref{eq:PDE2} calculated for two parameter sets PDE$_1$, PDE$_2$ plotted together with the CW solution of the DAE \eqref{eq:QGTI1}-\eqref{eq:QGTI3}. Although the bistability is lost for  PDE$_1$, it reappears for the  adjusted parameter values.
Parameters for PDE$_1$, DAE$_1$ are as in Figs.~\ref{fig:fig1}, \ref{fig:fig3}; for PDE$_2$ $u=0.52$, $\delta=0.044$ and $\gamma=10$.\label{fig:fig6}}
\end{figure}
After inserting all expansions into the NDDE~\eqref{eq:NDDE}, one collects orders of $\epsilon$ and solves them separately imposing the proper solvability conditions, see Supplementary Material for details. Finally, we find that the dispersive master PDE for the field $E$ on a fast timescale $\sigma$ containing the dynamics within one round-trip and an effective slow timescale $\theta$ that describes the dynamics over many round-trips, reads
\begin{equation}
 \partial_\theta E = \left(-i\,d_2\,\partial_\sigma^2+d_3\,\partial_\sigma^3\right)E+\left(\mathcal{L}+\mathcal{N}D\right)\,E+h\,Y_0\frac{\sqrt{1-\eta^2}}{1+i\delta}\,.\label{eq:PDE1}
\end{equation}
Here, the values of the coefficients $d_2$ and $d_3$ of the linear spatial
operator in Eq.~\eqref{eq:PDE1} are particularly instructive and read
\begin{equation}
 d_2=\frac{2\,\delta}{(1+\delta^2)^2}\,,\qquad d_3=\frac{2}{3}\frac{3\,\delta^2-1}{(1+\delta^2)^3}\,,
\end{equation}
whereas the coefficients $\mathcal{L}$ and $\mathcal{N}$ are
\begin{align}
    \mathcal{L}=&\,\eta-1-2\,\eta\,\arctan(\delta)\,i+2\,f\,\left(\frac{1-i\delta}{1+\delta^2}\right)^2\chi_0\,i,  \nonumber\\
    \mathcal{N}=&\,2\,f\,\frac{1-i\delta}{1+\delta^2}\,\chi^\prime_0\,i\,, \label{eq:PDE_coef}
\end{align}
Further, the equation for the carrier density $D$ takes the form
\begin{equation}
 \gamma^{-1}\,\partial_\sigma D=-D+\Im{\left(\chi_0+\chi^\prime_0D\right)}|E|^2\,.\label{eq:PDE2}
\end{equation}
One can see that the second order dispersion coefficient $d_2$ in Eq.~\eqref{eq:PDE1} cancels at resonance, i.e., at $\delta=0$, corresponding to the transition from anomalous to normal dispersion regime. Hence, at small $\delta$ values the dynamics is mostly controlled by TOD, which explains the oscillating tails of the TLSs observed in Fig.~\ref{fig:fig3}. Further, TOD vanishes for $\delta=\delta_c=\pm1/\sqrt{3}$.

Note that the linear  operator in Eq.~\eqref{eq:PDE1}  corresponding to the case of the linear empty cavity (i.e., $D=0$ (and $\chi=0$)) coincides with the linear part of the normal form PDE derived for the case of a quasi-conservative dispersive microcavity containing a Kerr medium coupled to a distant external mirror~\cite{SGJ-PRL-22}. There, a similar multiscale analysis performed in the long delay limit for small cavity losses and weak injection yielded the Lugiato-Lefever equation with TOD~\cite{LL-PRL-87,PGL-OL-14,PGG-PRA-17}. The latter can be understood as a dissipative version of the nonlinear Schr\"odinger equation in which the phase symmetry is broken by the presence of monochromatic injection with amplitude $Y_0$ and frequency offset $\delta$ with respect to the microcavity resonance. In the presence of the nonlinear element, the corresponding  contributions depending on $\chi_0$, $\chi'_{0}$, $f$ as well as $\delta$  appear in the coefficients $\mathcal{L}$, $\mathcal{N}$, cf. Eqs.~\eqref{eq:PDE_coef}.
\begin{figure}
\includegraphics[width=1\columnwidth]{fig7.pdf}
\caption{\label{fig:fig7} (a): Space-time evolution of a bright TLS obtained from the numerical simulations of the PDE model~\eqref{eq:PDE1}-\eqref{eq:PDE2}  for  $Y_0=0.6437$ (b),(c): Exemplary profiles taken at the last time-step showing the intensity  $|E|^2$ and $D$, respectively.}
\end{figure}

The CW solution $(E_s,\, D(E_s))$ of Eqs.~\eqref{eq:PDE1}-\eqref{eq:PDE2} can be found analytically by solving them for $Y_0 $ yielding
\begin{align}
    Y_0^2&=\frac{1+\delta^2}{h^2(1-\eta^2)}|\mathcal{L}+\mathcal{N}D(I_s)|^2\,I_s\,, \label{eq:CW1}\\
    D(I_s)&=\frac{\Im{(\chi_0)}\,I_s}{1-\Im{(\chi_0^\prime)}\,I_s}\,,\label{eq:CW2}
\end{align}
where $I_s=|E_s|^2$.
The resulting CW states \eqref{eq:CW1}-\eqref{eq:CW2}  calculated for two parameter sets labeled PDE$_1$ and PDE$_2$, respectively, are depicted together in Fig.~\ref{fig:fig6} with the CW solution DAE$_1$ of the DAE system ~\eqref{eq:QGTI1}-\eqref{eq:QGTI3} corresponding to the solution presented in Fig.~\ref{fig:fig3}. One can see that the PDE$_1$ solution for the same parameter set as the CW solution DAE$_1$ does not exhibit the bistability region. However, for a slightly different set of system parameter where the parameters that were scaled or expanded during the PDE~ \eqref{eq:PDE1}-\eqref{eq:PDE2} derivation are adjusted, the bistability region  is established  and is almost identical to the results obtained from the  DAE system ~\eqref{eq:QGTI1}-\eqref{eq:QGTI3}.
This discrepency can be ascribed to the fact that asumption of a small population inversion $D\ll 1$ is not always strictly verified during the emission of a TLS, see for instance Fig.~\ref{fig:fig7}~(c) and Fig.~\ref{fig:fig8}~(c).
\begin{figure}
\includegraphics[width=1\columnwidth]{fig8.pdf}
\caption{\label{fig:fig8} (a): Space-time evolution of a dark TLS obtained from the numerical simulations of the PDE model~\eqref{eq:PDE1}-\eqref{eq:PDE2}  for $Y_0=0.6491$. (b),(c): Exemplary profiles taken at the last time-step showing the intensitiy  $|E|^2$ and $D$, respectively.}
\end{figure}

In order to compare the solutions of the normal form  \eqref{eq:PDE1}-\eqref{eq:PDE2}  with
those obtained from the full DAE model~\eqref{eq:QGTI1}-\eqref{eq:QGTI3}, we solved the PDE master equations  \eqref{eq:PDE1}-\eqref{eq:PDE2} numerically for the parameter set PDE$_2$. The direct numerical simulations were performed using a split-step based pseudo-spectral method. The carrier equation~\eqref{eq:PDE2} was integrated using a dynamical boundary condition as in~ \cite{Hausen20}.

 To visualize the evolution of the resulting solution over many round-trips, a space-time representation in $(\sigma,\,\theta)$ plane is employed and the results are presented in Figs. ~\ref{fig:fig7},~\ref{fig:fig8}. Figure \ref{fig:fig7}~(a) shows the time evolution of a bright TLS.  The corresponding profiles at the end of the simulations are presented in panels (b) and (c) for the intensity $|E|^2$ and carrier density $D$, respectively. Here, one can see that  the domain walls lock at the second locking distance and the corresponding profile agrees with the bright TLS found in the DAE system, cf. Fig~\ref{fig:fig3}~(c).
An example of a  dark TLS is shown in Fig.~\ref{fig:fig8}~(a). The initial condition evolves into a dark TLS with strong asymmetric oscillatory tail as seen in an exemplary profiles in (b) and (c) for $E$ and $D$, respectively. Again, this solution is in good agreement with the dark TLS found in the original DAE model, see Fig~\ref{fig:fig3}~(b).
Note that both dark TLSs from the PDE and DAE possess more oscillatory tails as the bright TLSs.  This situation is similar to the one observed for the Lugiato-Lefever equation subjected to TOD in the normal dispersion regime~\cite{PGG-PRA-17}. There, TOD allows for stable dark and bright TLSs to coexist, while the shape of the switching waves connecting the top and bottom CW solutions is modified and bright TLSs are stabilized due to the generation of additional oscillations in the switching-wave profiles.

\section{Square-Waves}
Finally, we consider the same system given by DAE~\eqref{eq:QGTI1}-\eqref{eq:QGTI3}  but operated out of the bistability regime such that the wavelength of the optical injection is set in-between two modes of the external cavity. This situation corresponds to the case of antiresonant feedback, i.e. $\varphi=\pi$.

Note that the formation of SWs was reported in various optical and opto-electronic systems subjected to time-delayed feedback~\cite{PJC-PRE-09,JAH-PRL-15,MGJ-PRA-07,GES-OL-06,LLZ-OL-16}.  There, they typically appear in the long delay limit  via a supercritical Andronov-Hopf (AH) bifurcation as modulated oscillations with a period close to twice the value of the time-delay $\tau$. The SWs then evolve into sharp transition layers connecting plateaus of duration $\sim\tau$. Recently, the formation of SWs in a monomode micro-cavity containing a Kerr medium coupled to a long external feedback cavity under CW  injection in the normal dispersion regime was demonstrated~\cite{KSG-OL-22}.  There,  SWs also appear via a homoclinic snaking scenario, leading to the formation of complex-shaped multistable SW solutions.
 \begin{figure}
	\includegraphics[width=1\columnwidth]{fig9.pdf}
	\caption{SWs emerging in a supercritical AH bifurcation at $H_1$ and a subcritical AH bifurcation at $H_2$. Such a branch will be referred to as a subcritical SW branch of the QGTI. b) shows a stable SW with $50\%$ duty cycle, while c) shows a solution with three plateaus at the Maxwell line.  $(\delta,h,\eta,\varphi,\gamma,f,u,b,\tau)=(-0.01,2,0.98,\pi,0.1,1,0.65,1000,1000)$. \label{fig:fig9}}
	\end{figure}
To study the formation and the properties of the SWs in the DAE model~\eqref{eq:QGTI1}-\eqref{eq:QGTI3},  we perform the path-continuation as a function of the injection $Y_0$,  setting the feedback phase $\varphi=\pi$. Our results are shown in Fig.~\ref{fig:fig9}. In Fig.~\ref{fig:fig9}~(a), for increasing $Y_0$, a branch of stable SW solutions (dark red) emerges from the stable CW state (black) at a supercritical AH bifurcation at the point $H_1$. Here, a small-amplitude  periodic solution with a period $\sim 2\tau$ close to the bifurcation point transforms into a well-developed SW with sharp transition layers  when $Y_0$ is increased, see Fig.~\ref{fig:fig9}~(b), where the exemplary profile is shown (cf. the blue point in the panel (a)). However, for further increased $Y_0$, the SW branch turns in a saddle-node bifurcation (cf. the inset in the panel (a)) and experiences a sequence of further saddle-node bifurcations  resulting in a collapsed snaking structure. In particular, close to  the vertical Maxwell line where the plateaus shrink in favor of the CW value, it leads to a three-plateau solution as shown in panel (c). Finally, the SW branch merges with the CW state at $H_2$ in a subcritical AH bifurcation. Remarkably, the snaking scenario is similar to one reported in ~\cite{KSG-OL-22}, where, however, it was observed for low injection values, whereas in the present case the subcritical AH bifurcation is created at high $Y_0$.

By varying $\eta$ to lower values one can observe a transition of the AH bifurcation at $H_2$ from subcritical to supercritical, which leads to a disappearance of the snaking structure. The transition seems to be located close to the SWs existing threshold in $\eta$. A resulting branch can be seen in Fig.~ \ref{fig:fig10}~(a). While panel (b) demonstrates an exemplary profile of a  fully developed SW with sharp transition layers, panel (c) shows a smoother periodic solution close to the (now) supercritical AH bifurcation at $H_2$.

  \begin{figure}
	\includegraphics[width=1\columnwidth]{fig10.pdf}
	\caption{An exemplary branch for the supercritical regime, close to the transition from sub- to supercritical SW solutions. The two profiles in b) and c) show a fully developed as well as a just emerging SW, respectively. The parameters are $(\delta,h,\eta,\varphi,\gamma,f,u,b,\tau)=(-0.01,2,0.86,\pi,0.1,1,0.65,1000,1000)$ \label{fig:fig10}.}
	\end{figure}

\section*{Conclusion}

In  conclusion, in this paper  the dynamics of temporal localized states and square-waves in semiconductor micro-resonators with strong time-delayed feedback was analyzed. We considered a semiconductor quantum-well as the nonlinear element instead of an ideal, instantaneous Kerr response. This allowed us to study the influence of the band-gap detuning with respect to the injection field and the population inversion time scale onto the dynamics. A first principle model for the optical response was used in combination with time delayed algebraic constraints that correspond to the boundary conditions linking the optical field in the various parts of the compound cavity system.

Using a combination of analytical, numerical and path-continuation methods, the existence of  dark and bright temporal localized states  appearing via the locking of domain walls that connect the high and low bistable continuous wave  intensity levels was demonstrated in the normal dispersion regime and for a feedback phase of zero. Further, operated in the antiresonant feedback regime, square-waves solutions, possessing a period of approximately twice the round-trip in the external cavity were obtained. Depending on the feedback strength that is connected with the reflectivity of the external mirror, the  latter can also appear through a homoclinic snaking scenario. Finally, to clarify the role of the second and third order dispersion in the formation of the localized states, a multiple time-scale analysis in the good cavity limit was performed. The resulting dispersive normal form partial differential equation shows a good qualitative agreement with the original time-delayed model.

More generally, the results presented here may serve as a guide for further analysis into the dynamics of  temporal localized states in semiconductor micro-resonators subjected to the strong time-delayed feedback in the long cavity regime. Our analysis clarified that, even in the good cavity limit, when the carrier lifetime becomes much smaller than the micro-cavity photon lifetime, both branches of bright and dark TLSs disappear, even if CW optical bistability is not affected. This indicates that the most favorable situation would certainly consider a quantum-well medium optimized to act as a saturable absorber with a fast carrier recovery time in the picosecond range instead of, e.g., a micro-cavity containing a gain medium engineered to favor a lasing action and for which the carrier lifetime would be in the nanosecond range.

\section*{Supplementary Material}
See supplementary material for the animation showing the solution profiles along the branches shown in Fig. \ref{fig:fig3} in Sec.~\ref{sec:DAE} and for details on the derivation of the normal form partial differential equation presented in Sec.~\ref{sec:PDE}.
\section*{Author Declarations}
The authors have no conflicts to disclose.
\section*{Acknowledgments}
T.S. thanks the foundation “Studienstiftung des deutschen Volkes” for financial support,  E.K, J. J. and S.G. acknowledge the financial support of the projects KEFIR, Ministerio de Economía y Competitividad (PID2021-128910NB-100 AEI/FEDER UE) and KOGIT, Agence Nationale de la Recherche (ANR-22-CE92-0009)), Deutsche Forschungsgemeinschaft (DFG) via Grant Nr. 505936983.
We thank A. Giacomotti and B. Garbin for useful discussions.
\section*{Data Availability}
The data that support the findings of this study are available from the corresponding author upon reasonable request.

\section*{References}

\section*{Supplementary material: Normal Form derivation}

In this supplementary material we demonstrate how a partial differential equation (PDE) normal form model can be derived from our first principle Delay-Algebraic Equation (DAE) for a micro-cavity operating in the Gires-Tournois regime containing a quantum well and coupled to an long external cavity under CW injection using multiple time scaling. The resulting PDE system can help to get a better understanding on how the parameters influence the dynamics of temporal localized states in the underlying system. The derivation will be done in a good cavity limit, i.e. $\eta\rightarrow 1$ and $h=2$.

\section{Normal Form Settings}
We start from the full DAE system which reads
\begin{align}
    \dot{E}(t)&=\left[if\chi -1 -i\delta\right]E(t)+hY(t),\label{eq:QGTI_1}\\
    \dot{D}(t)&=\gamma\left[-D(t)+\Im(\chi)|E|^2\right],\label{eq:QGTI_2}\\
    Y(t)&=\eta e^{i\varphi}\left[E(t-\tau)-Y(t-\tau)\right]+Y_0\sqrt{1-\eta^2}, \label{eq:QGTI_3}
\end{align}
with
\begin{align}
    \chi(u,D)=\frac{1}{\pi}\left(\ln(u+i-b)-2\ln(u+i-D)+\ln(u+i)-\ln(b)\right). \label{eq:QGTI_chi}
\end{align}
Here, we only consider the case $\varphi=0$. Going from a DAE to a neutral delay differential equation (NDDE) allows to eliminate $Y(t)$. This can be achieved by adding Eq. \eqref{eq:QGTI_1} at time $t$ onto itself at time $t-\tau$, multiplied by $\eta$. One can then recognize and insert Eq. \eqref{eq:QGTI_3}, which eliminates $Y(t)$ and leads to
\begin{align}
    \left(\dot{E}(t)+\eta\dot{E}(t-\tau)\right)=&\left[if\chi(D(t))-1-i\delta\right]E(t)+\eta\left[if\chi(D(t-1))-1-i\delta\right]E(t-\tau)\nonumber\\&+h\eta E(t-\tau)+hY_0\sqrt(1-\eta^2).
\end{align}
\subsection{Time Scaling}
First, the system is rescaled with $t\rightarrow\frac{t}{\tau}$ and $\frac{1}{\tau}=\epsilon$, yielding
\begin{align}
    \epsilon\left(\dot{E}(t)+\eta\dot{E}(t-1)\right)=&\left[if\chi(D(t))-1-i\delta\right]E(t)+\eta\left[if\chi(D(t-1))-1-i\delta\right]E(t-1)\nonumber\\
    &+h\eta E(t-1)+hY_0\sqrt(1-\eta^2)\label{eq:rescaled_QGTI1_3},\\
   \epsilon\dot{D}=& \gamma\left[-D(t)+\Im(\chi(D))|E|^2\right].
\end{align}
It now has a period of $T \approx 1$ defining the fastest timescale $t_0$. To avoid drift, the period $T$ or equivalently the frequency $\omega$ can be expanded allowing for higher order corrections. This leads to a stretch in $t_0$:
\begin{align}
    t_0&=\omega t,\\
    t_1&=0,\\
    t_2&=\epsilon^2t,\\
    t_3&=\epsilon^3t.
\end{align}
\subsection{Parameter Expansions}
In the following, the parameters will be expanded based on the knowledge about the system as well as some practical thoughts:
\begin{align}
    Y_0&=\epsilon\gamma_1+\epsilon^2\gamma_2+\epsilon^3\gamma_3,\\
    \eta&=1+\eta^2\eta_2.
\end{align}
As mentioned before, the normal form is derived in a good cavity limit by choosing $\eta\rightarrow1$ and $h=2$ eliminating $h$ from the equation. Furthermore, choosing $Y_0$ being $\mathcal{O}(\epsilon)$ combined with the choice of $\eta$  leads to an injection that comes in at $\mathcal{O}(\epsilon^2)$. The detuning $\delta$ will not be expanded, nevertheless $\delta$ will be assumed to be small at some point of the derivation. The parameters $u$ and $b$ only appear within $\chi$, hence they will not be considered for the moment.
Finally, the major parameter to be discussed is $f$ determining the strength of the nonlinearity.
\subsection{Nonlinearity}
In general the nonlinearity should enter at $\mathcal{O}(\epsilon^3)$ limiting nonlinear effects. The most obvious path is the choice $f=\mathcal{O}(\epsilon^3)$. In that case however, it appears to be hard to access the nonlinearity itself from the argument of $\chi$. Another way is to choose $f$ in a lower order  such as $f\rightarrow \epsilon^2 f$ and rescale the field, as the following example will demonstrate. Consider $E\rightarrow\epsilon^\frac{1}{2} \mathcal{E}$, $Y\rightarrow\epsilon^\frac{1}{2} \mathcal{Y}$ and $Y_0\rightarrow \epsilon^\frac{1}{2} \mathcal{Y}_0$. This leads to:
\begin{align}
    \epsilon \dot{D}(t)&=\gamma\left[-D(t)+\epsilon\Im(\chi(D))|\mathcal{E}|^2\right].
\end{align}
Now $D\rightarrow \epsilon \mathcal{D}$ is an obvious consequence to be able to divide each term by $\epsilon$ yielding
\begin{align}
     \epsilon\dot{\mathcal{D}}(t)&=\gamma\left[-\mathcal{D}(t)+\Im(\chi(\epsilon \mathcal{D}))|\mathcal{E}|^2\right],
\end{align}
as well as $\chi(\epsilon \mathcal{D})$ in general. For Eq. \eqref{eq:rescaled_QGTI1_3} one finds
\begin{align}
    \epsilon^\frac{1}{2}\epsilon\left(\dot{\mathcal{E}}(t)+\eta\dot{\mathcal{E}}(t-1)\right)=&\epsilon^\frac{1}{2}\left[i\epsilon^2f\chi(\epsilon \mathcal{D})-1-i\delta\right]\mathcal{E}(t)\nonumber\\
    &+\epsilon^\frac{1}{2}\eta\left[i\epsilon^2f\chi(\epsilon \mathcal{D}(t-1))-1-i\delta\right]\mathcal{E}(t-1)\nonumber\\
    &+\epsilon^\frac{1}{2}h\eta \mathcal{E}(t-1)+\epsilon^\frac{1}{2} h\mathcal{Y}_0\sqrt(1-\eta^2). \label{eq:insert1}
\end{align}
Here $\epsilon^\frac{1}{2}$ can be taken away from each term to reproduce the old equation only differing in the argument of $\chi$. To write the susceptibility as a polynomial in $\epsilon$, a Taylor expansion can be included for small $\epsilon$, leading to:
\begin{align}
    \chi(\epsilon \mathcal{D})=&\chi(0)+\epsilon \chi'(0)\mathcal{D} + \mathcal{O}(\epsilon^2).
\end{align}
Note that the higher order terms will not occur in the final equation if only a total order of $\mathcal{O}(\epsilon^3)$ is considered. The terms $\chi(0)$ and $\chi'(0)$ are the complex terms
\begin{align}
   \chi(0)&=\frac{1}{\pi}\ln\left(\frac{u+i-b}{(u+i)b}\right)=\frac{1}{\pi}\ln\left(\frac{1}{b}-\frac{1}{u+i}\right)\approx\frac{1}{\pi}\ln\left(\frac{-1}{u+i}\right)=-\frac{1}{\pi}\ln(-u-i),\\
   \chi'(0)&=\frac{1}{\pi}\frac{2}{u+i-0}=\frac{1}{\pi}\frac{2(u-i)}{u^2+1},
\end{align}
which do not depend on the fields. Here, $b\gg 1$ has been used. This approximation is optional and can be omitted easily, if needed, since the susceptibility will be inserted in the final equation. The imaginary part of $\chi(0)$ can be isolated using the Arg function:
\begin{align}
    -\ln(-u-i)=&-\ln|u+i|-i\text{Arg}(-u-i)\\
    =&-\ln\sqrt{u^2+1}+i\left[\frac{\pi}{2}+\arctan(u)\right].
\end{align}
This leads to
\begin{align}
    \chi(0)&=-\frac{1}{\pi}\ln\sqrt{u^2+1}+i\left[\frac{1}{2}+\frac{1}{\pi}\arctan(u)\right]=R_0+iI_0,\\
    \chi'(0)&=\frac{2}{\pi}\frac{u-i}{u^2+1}=R_1+iI_1,
\end{align}
which yields an approximated susceptibility
\begin{align}
    \chi(\epsilon \mathcal{D})=R_0+iI_0+\epsilon(R_1+iI_1)\mathcal{D}+\mathcal{O}(\epsilon^2),
\end{align}
with the nonlinear term in $\mathcal{O}(\epsilon^3)$. Inserting these results one finds for $\dot{\mathcal{D}}$
\begin{align}
   \epsilon  \dot{\mathcal{D}}(t)&=\gamma\left[-\mathcal{D}(t)+\Im(\chi(\epsilon \mathcal{D}))|\mathcal{E}|^2\right]\\
     &=\gamma\left[-\mathcal{D}(t)+(I_0+\epsilon I_1 \mathcal{D}(t))|\mathcal{E}(t)|^2\right]. \label{eq:insert2}
\end{align}

\section{Derivation}
\subsection{Time Derivative and Drift Operator}
Now the corresponding time derivative operator is:
\begin{align}
    \partial_t=\omega\partial_0+\epsilon^2\partial_2+\epsilon^3\partial_3.
\end{align}
To equalize the drift the expression $\omega=1+\epsilon\omega_1+\epsilon^2\omega_2+\epsilon^2\omega_3$ can be used which leads to:
\begin{align}
    \partial_t&=\partial_0+\epsilon\omega_1\partial_0+\epsilon^2(\omega_2\partial_+\partial_2)+\epsilon^3(\omega_3\partial_0+\partial_3)+\mathcal{O}(\epsilon^4)\\
    &=\mathcal{T}_0+\epsilon\mathcal{T}_1+\epsilon^2\mathcal{T}_2+\epsilon^3\mathcal{T}_3.
\end{align}
The field has to be written with respect to the new time scale. Thus, $\mathcal{E}(t)=\mathcal{E}(t_0,t_2,t_3)$, but also $\mathcal{E}(t-1)=\mathcal{E}(t_0-\omega,t_2-\epsilon^2,t_3-\epsilon^3)$ has to be used. In order to sort the system by orders of $\epsilon$ an expansion for small $\epsilon$ can be used to rewrite the time-delayed field
\begin{align}
    \mathcal{E}(t_0-\omega,t_2-\epsilon^2,t_3-\epsilon^3)=&\mathcal{E}(t_0-1-\epsilon\omega_1-\epsilon^2\omega_2-\epsilon^3\omega_3,t_2-\epsilon^2,t_3-\epsilon^3)\\
    =&\mathcal{E}(t_0-1,t_2,t_3)+\epsilon(-\omega_1\partial_0)\mathcal{E}(t_0-1,t_2,t_3)\nonumber\\
    &+\epsilon^2\left(\frac{\omega_1^2\partial_0^2}{2}-\omega_2\partial_0-\partial_2\right)\mathcal{E}(t_0-1,t_2,t_3)\\\nonumber&+\epsilon^3\left(\frac{-\omega_1^3\partial_0^3}{6}+\omega_1\omega_2\partial_0^2+\omega_1\partial_0\partial_2-\omega_3\partial_0-\partial_3\right)\mathcal{E}(t_0-1,t_2,t_3)
\end{align}
The found prefactors correspond to the shift operator $\mathcal{S}$ fulfilling
\begin{align}
    \mathcal{E}(t_0-\omega,t_2-\epsilon^2,t_3-\epsilon^3)&=\mathcal{S}\mathcal{E}(t_0-1)\\
    &=\left(1+\epsilon\mathcal{S}_1+\epsilon^2\mathcal{S}_2+\epsilon^3\mathcal{S}_3\right)\mathcal{E}(t_0-1),
\end{align}
with the reduced notation $\mathcal{E}(t_0-1,t_2,t_3)=\mathcal{E}(t_0-1)$.\\
Using the operators introduced above all definitions can be inserted into Eq. \eqref{eq:insert1} and Eq. \eqref{eq:insert2}. Additionally, the fields are expanded in a power series as $\mathcal{E}=\mathcal{E}_0+\epsilon\mathcal{E}_1+\epsilon^2\mathcal{E}_2+\epsilon^3\mathcal{E}_3+\mathcal{O}(\epsilon^4)$. Afterwards the resulting expressions will be sorted in orders of $\epsilon$ which will be considered separately.

\subsection{Solving by Orders}
Starting with the lowest order $\mathcal{O}(\epsilon^0)$ one obtains
\begin{align}
\mathcal{E}_0(t_0)-\mathcal{E}_0(t_0-1)=0,
\end{align}
describing a simple periodic solution. \\
Using this, from $\mathcal{O}(\epsilon^1)$
\begin{align}
\mathcal{E}_1(t_0)-\mathcal{E}_1(t_0-1)&=\left(-\omega_1-\frac{2}{\delta^2+1}\right)\partial_0\mathcal{E}_0(t_0),
\end{align}
is found.
This equation describes a periodic solution with a periodicity disturbed by a drift term elucidating the meaning of the expansion in $\omega$ as compensation for potential drift. Without the first order term $\omega_1\partial_0\mathcal{E}_0(t_0)$ a nonvanishing drift of $\mathcal{E}_0$ on the time scale $t_0$ is observed. The drift is compensated through the result for the small frequency correction $\epsilon\omega_1$ by choosing
\begin{align}
    \omega_1=\frac{-2}{1+\delta^2}.
\end{align}
For $\mathcal{O}(\epsilon^2)$, the lower orders will be used once again leading to $\omega_2=\omega_1^3$. Furthermore, the solvability condition
\begin{align}
    \partial_2\mathcal{E}_0=\left[-i\delta\frac{\omega_1^2}{2}\partial_0^2-L_2\right]\mathcal{E}_0+\Gamma_2,
\end{align}
is found with
\begin{align}
    L_2&=\frac{2(fI_0-ifR_0)}{(1+i\delta)^2}-(\eta_2+i\theta_2),\\
    \Gamma_2&=\gamma_1\frac{2\sqrt{-2\eta_2}}{1+i\delta}.
\end{align}
The operator $L_2$ contains the first term from the expansion of the susceptibility, but the nonlinearity itself does not enter until $\mathcal{O}(\epsilon^3)$.
Here, $\omega_3=\omega_1^3$ is found.
The final solvability condition for the third order containing all simplifications, has the form
\begin{align}
    \partial_3\mathcal{E}_0+\partial_2\mathcal{E}_1=&-i\frac{3\delta}{2}\omega_1^3\partial_0^2\mathcal{E}_0+\frac{1-3\delta^2}{12}\omega_1^3\partial_0^3\mathcal{E}_0\nonumber\\
    &+[\mathcal{N}\mathcal{D}-\omega_1L_2]\mathcal{E}_0+\left[-i\delta\frac{\omega_1^2}{2}\partial_0^2-L_2\right]\mathcal{E}_1+\Gamma_3+\omega_1\Gamma_2,
\end{align}
with
\begin{align}
    L_2&=\frac{2(fI_0-ifR_0)}{(1+i\delta)^2}-(\eta_2+i\theta_2),\\
    \Gamma_2&=\gamma_1\frac{2\sqrt{-2\eta_2}}{1+i\delta},\\
    \Gamma_3&=\gamma_2\frac{2\sqrt{-2\eta_2}}{1+i\delta},\\
    \mathcal{N}&=\frac{2f(iR_1-I_1)}{1+i\delta}.
\end{align}
Equation \eqref{eq:QGTI_2} of the $\mathcal{D}$ field that has been rewritten to
\begin{align}
   \epsilon  \dot{\mathcal{D}}(t)=\gamma\left[-\mathcal{D}(t)+(I_0+\epsilon I_1 \mathcal{D}(t))|\mathcal{E}(t)|^2\right],
\end{align}
has to be discussed as well. Inserting the time derivative operator as well as the nonlinearity leads to
\begin{align}
    \left(\epsilon\partial_0+\epsilon^2\omega_1\partial_0+\epsilon^3(\omega_2\partial_0+\partial_2)\right) \dot{D}(t)
     &=\gamma\left[-D(t)+(I_0+\epsilon I_1 D(t))|E(t)|^2\right].
\end{align}
    \section{Assembling the Full Partial Differential Equation System}
Finally the PDE shall be deduced from the previous results. Therefor, an effective time scale
\begin{align}
    \partial_\theta=\epsilon^2\partial_2+\epsilon^3\partial_3
\end{align}
is introduced. The results are now included by using $\mathcal{E}=\mathcal{E}_0+\epsilon \mathcal{E}_1$, where higher orders of the field fall into $\mathcal{O}(\epsilon^4)$ due to the definition of the time scale. This leads to:
\begin{align}
    \partial_\theta \mathcal{E}
    =&\epsilon^2\left(\left[-i\delta\frac{\omega_1^2}{2}\partial_0^2-L_2\right]\mathcal{E}+(\gamma_1+\epsilon\gamma_2)\frac{2\sqrt{-2\eta_2}}{1+i\delta}\right)\nonumber\\
    &+\epsilon^3\left(\frac{1-3\delta^2}{12}\omega_1^3\partial_0^3\mathcal{E}+[\mathcal{N}\mathcal{D}-\omega_1L_2]\mathcal{E}+\omega_1\gamma_1\frac{2\sqrt{-2\eta_2}}{1+i\delta}\right).
\end{align}
Next, space can be rescaled as $\epsilon|\omega_1|\partial_0\rightarrow\partial_\sigma$, where the absolute value leads to a switch in sign of the third order dispersion:
\begin{align}
    \partial_\theta \mathcal{E}=&\left[\frac{-i\delta}{2}\partial_\sigma^2-\epsilon^2(-\eta_2+2W_1-i\theta_2)\right]\mathcal{E}+h\mathcal{Y}_0\frac{\sqrt{1-\eta^2}}{1+i\delta}-\frac{1-3\delta^2}{12}\partial_\sigma^3\mathcal{E}\nonumber\\
    &+\epsilon^3\left[\frac{2f(iR_1-I_1)}{1+i\delta}\mathcal{D}-\omega_1(-\eta_2+2W_1-i\theta_2)\right]\mathcal{E}+\epsilon\omega_1h\mathcal{Y}_0\frac{\sqrt{1-\eta^2}}{1+i\delta}+\mathcal{O}(\epsilon^4).
\end{align}
Applying the parameter expansions yields $\epsilon^2(\eta_2+i\theta)+\mathcal{O}(\epsilon^4)=\eta e^{i\theta}-1$. Also, $\epsilon^2f \rightarrow f$ is reset and small terms containing $\epsilon\omega_1$ are removed, yielding
\begin{align}
    \partial_\theta \mathcal{E}=&\left[\frac{-i\delta}{2}\partial_\sigma^2-2f\frac{I_0-iR_0}{(1+i\delta)^2}-1+\eta e^{i\theta}\right]\mathcal{E}+h\mathcal{Y}_0\frac{\sqrt{1-\eta^2}}{1+i\delta}\nonumber\\
    &-\frac{1-3\delta^2}{12}\partial_\sigma^3\mathcal{E}+\epsilon\frac{2f(iR_1-I_1)}{1+i\delta}\mathcal{D}\mathcal{E}+\mathcal{O}(\epsilon^4).
\end{align}
One can progress by multiplying both sides of the equation by $\epsilon^\frac{1}{2}$ and inverting the field's scaling with $\mathcal{E}\epsilon^\frac{1}{2}\rightarrow E$, $\mathcal{Y}_0\epsilon^\frac{1}{2}\rightarrow Y_0$ and $\mathcal{D}\epsilon\rightarrow D$, which results in
\begin{align}
    \partial_\theta E=&\frac{-i\delta}{2}\partial_\sigma^2E-\frac{1-3\delta^2}{12}\partial_\sigma^3E+\frac{2f(iR_1-I_1)}{1+i\delta}DE\nonumber\\&
    +\left[-2 f\frac{I_0-iR_0}{(1+i\delta)^2}-1+\eta e^{-i2\arctan(\delta)}\right]E
    + hY_0\frac{\sqrt{1-\eta^2}}{1+i\delta}.
\end{align}
Here, the original fields from Eqs. \eqref{eq:QGTI_1} - \eqref{eq:QGTI_chi} are recovered. This is necessary to compare the derived PDE with the underlying DAE system.\\
After expanding the exponential function for small values of $\delta$ only the coefficients $R_0$, $I_0$, $R_1$, and $I_1$ remain. Hence, the approximation for the susceptibility shall be finally inserted leading to
\begin{align}
    \partial_\theta E=&\frac{-i\delta}{2}\partial_\sigma^2E-\frac{1-3\delta^2}{12}\partial_\sigma^3E+\frac{4f}{\pi}\frac{(1+iu)(1-i\delta)}{(1+u^2)(1+\delta^2)}DE\nonumber\\
    &+\left[-\frac{2f}{\pi}(1-i\delta)^2\frac{i\ln\sqrt{u^2+1}+\frac{\pi}{2}+\arctan(u)}{\left(1+\delta^2\right)^2}-1+\eta-i2\eta\arctan(\delta)\right]E\nonumber\\
    &+ hY_0\frac{\sqrt{1-\eta^2}}{1+i\delta},
\end{align}
which is this section's final PDE for $E$. The second equation
\begin{align}
    (\epsilon\partial_0+\epsilon^2\omega_1\partial_0+\epsilon^3(\omega_2\partial_0+\partial_2)) \mathcal{D}(t)
     &=\gamma\left[-\mathcal{D}(t)+(I_0+\epsilon I_1 \mathcal{D}(t))|\mathcal{E}(t)|^2\right]
\end{align}
will be simplified using $\epsilon|\omega_1|\partial_0=\partial_\sigma$ as well as $\epsilon^2\partial_2+\omega^3\partial_3=\partial_\theta$. However, one needs to keep in mind that the fields involved still have to be rescaled. Thus, $\mathcal{E}\epsilon^\frac{1}{2}\rightarrow E$ and $\mathcal{D}\epsilon\rightarrow D$ are applied. This leads to
\begin{align}
    \left(\frac{1}{-\omega_1}\partial_\sigma-\epsilon(\partial_\sigma+\partial_\theta)+\epsilon^3\omega_2\partial_0+\omega^4\omega_3\partial_0\right) \frac{1}{\epsilon} D
     &=\gamma\left[\frac{-1}{\epsilon} D+(I_0+ I_1\frac{\epsilon}{\epsilon} D)\frac{1}{\epsilon}|E|^2\right],
 \end{align}
 which is equivalent to
 \begin{align}\left(\frac{1}{-\omega_1}\partial_\sigma-\epsilon\partial_\sigma+\epsilon\partial_\theta-\epsilon^2\omega_1\partial_\sigma-\epsilon^3\omega_1^2\partial_\sigma\right) \frac{1}{\epsilon}D&=
     \frac{\gamma}{\epsilon}\left[-D+(I_0+I_1D)|E|^2\right].
\end{align}
Here, the relations $\omega_2=\omega_1^2$ and $\omega_3=\omega_1^3$ have been used to rewrite all derivative operators with the new scales $\sigma$ and $\theta$.
Now each side is multiplied by $\epsilon$. In a first attempt the lowest order can be used to describe the dynamics of $D$ as
\begin{align}
          \partial_\sigma D
     =-\omega_1\gamma\left[-D+(I_0+I_1D)|E|^2\right].
\end{align}
However, interpreting $\sigma$ as a quasi-spatial variable and $\theta$ as the time variable this is not a differential equation in time $\sigma$. Indeed this section leads to a PDE in $E(t)$ with a constraint for the field $D(t)$. The resulting system is

\begin{align}
    \partial_\theta E=&\frac{-i\delta}{2}\partial_\sigma^2E-\frac{1-3\delta^2}{12}\partial_\sigma^3E+\frac{4f}{\pi}\frac{(1+iu)(1-i\delta)}{(1+u^2)(1+\delta^2)}DE\nonumber\\
    &+\left[-\frac{2f}{\pi}(1-i\delta)^2\frac{i\ln\sqrt{u^2+1}+\frac{\pi}{2}+\arctan(u)}{\left(1+\delta^2\right)^2}-1+\eta-i2\eta\arctan(\delta)\right]E\nonumber\\
    &+ hY_0\frac{\sqrt{1-\eta^2}}{1+i\delta},\\
  0=&-\partial_\sigma D
     +\frac{\gamma}{1+\delta^2}\left[-2D+\left(1+\frac{2}{\pi}\arctan(u)-\frac{4}{\pi}\frac{1}{u^2+1}D\right)|E|^2\right].
\end{align}
The equations are finally simplified to
\begin{align}
    \partial_\theta E &= \left(-i \tilde{d}_2\partial_\sigma^2+\tilde{d}_3\partial_\sigma^3+\mathcal{L}\right)E+\mathcal{N}DE+hY_0\frac{\sqrt{1-\eta^2}}{1+i\delta},\label{eq:PDE1}\\
    0&=\left(\zeta_3\partial_\sigma-2-\zeta_2|E|^2\right)D+\zeta_1|E|^2,\label{eq:PDE2}
\end{align}
with
\begin{align}
    \tilde{d}_2&=\frac{\delta}{2},\\
    \tilde{d}_3&=\frac{3\delta^2-1}{12},\\
    \mathcal{L}&=-\frac{2f}{\pi}(1-i\delta)^2\frac{i\ln\sqrt{u^2+1}+\frac{\pi}{2}+\arctan(u)}{\left(1+\delta^2\right)^2}-1+\eta-i2\eta\arctan(\delta),  \\
    \mathcal{N}&=\frac{4f}{\pi}\frac{(1+iu)(1-i\delta)}{(1+u^2)(1+\delta^2)},\\
    \zeta_1&=\left(1+\frac{2}{\pi}\arctan(u)\right),\\
    \zeta_2&=\frac{4}{\pi}\frac{1}{u^2+1},\\
    \zeta_3&=-\frac{1+\delta^2}{\gamma}.
\end{align}

In the letter we used an scaled $\sigma$ scale to allow for a better comparison with previous results. Using $\sigma\rightarrow-\omega_1\sigma=\frac{2}{1+\delta^2}\sigma$ leads to
\begin{equation}
 \partial_\theta E = \left(-i\,d_2\,\partial_\sigma^2+d_3\,\partial_\sigma^3\right)E+\left(\mathcal{L}+\mathcal{N}D\right)\,E+h\,Y_0\frac{\sqrt{1-\eta^2}}{1+i\delta}\,.\label{eq:PDE1}
\end{equation}
Here, the values of the coefficients $d_2$ and $d_3$ of the linear spatial
operator in Eq.~\eqref{eq:PDE1} read
\begin{equation}
 d_2=\frac{2\,\delta}{(1+\delta^2)^2}\,,\qquad d_3=\frac{2}{3}\frac{3\,\delta^2-1}{(1+\delta^2)^3}\,,
\end{equation}
Additionally rewriting the equations as a function of $\chi(0)$ and $\chi^\prime(0)$ leads to the coefficients $\mathcal{L}$ and $\mathcal{N}$
\begin{align}
    \mathcal{L}=&\,\eta-1-2\,\eta\,\arctan(\delta)\,i+2\,f\,\left(\frac{1-i\delta}{1+\delta^2}\right)^2\chi_0\,i,  \nonumber\\
    \mathcal{N}=&\,2\,f\,\frac{1-i\delta}{1+\delta^2}\,\chi^\prime_0\,i\,, \label{eq:PDE_coef}
\end{align}
In that case, the equation for the carrier density $D$ takes the form
\begin{equation}
 \gamma^{-1}\,\partial_\sigma D=-D+\Im{\left(\chi_0+\chi^\prime_0D\right)}|E|^2\,.\label{eq:PDE2},
\end{equation}
which is the form presented in the letter.
\end{document}

\section{Literatur}}
\usepackage[top=3cm,bottom=3cm]{geometry}
\setlength\parindent{0pt}
\usepackage{csquotes}

\addto\extrasngerman{
	\let\subsectionautorefname\sectionautorefname
	\let\subsubsectionautorefname\sectionautorefname
}

\usepackage{chngcntr}
\counterwithin{figure}{section}
\counterwithin{table}{section}
\counterwithin{equation}{section}

\ifdraft{
	\usepackage{refcheck}
	\newcommand{\refcheckize}[1]{%
	  \expandafter\let\csname @@\string#1\endcsname#1%
	  \expandafter\DeclareRobustCommand\csname relax\string#1\endcsname[1]{%
	    \csname @@\string#1\endcsname{##1}\wrtusdrf{##1}}%
	  \expandafter\let\expandafter#1\csname relax\string#1\endcsname
	}
	\refcheckize{\autoref}
	\usepackage[german,colorinlistoftodos]{todonotes}
	\usepackage{draftwatermark}
	\SetWatermarkText{ENTFWURF}
	\SetWatermarkScale{5}
	\SetWatermarkColor[rgb]{0.7,0,0}
	\SetWatermarkLightness{0.9}
	\newcommand{\todoDiscuss}[1]{\todo[color=green!40]{#1}}
	\newcommand{\todoMissing}[1]{\todo[inline,color=red!40]{Fehlt noch: #1}}
	\newcommand{\todoNote}[1]{\todo[color=orange!40]{#1}}
	\newcommand{\todoMissingFigure}[1]{\missingfigure{#1}}
}{
	\newcommand{\todoDiscuss}[1]{\errmessage{Noch zu Diskutieren: #1}}
	\newcommand{\todoMissing}[1]{\errmessage{Fehlt noch: #1}}
	\newcommand{\todoNote}[1]{\errmessage{Noch zu tun: #1}}
	\newcommand{\todoMissingFigure}[1]{\errmessage{Bild Fehlt: #1}}
}

    \usepackage{siunitx}
    \numberwithin{equation}{section}
    \usepackage{braket}
    \usepackage{dsfont}
    \usepackage{accents}
    \usepackage{mathtools}