11institutetext: 1\. Munster Technological University, Cork, Ireland
2\. University College Cork, Cork, Ireland
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>
# Inter and Intra Signal Variance in Feature Extraction and Classification of
Affective State††thanks: This publication has emanated from research supported
in part by a Grant from Science Foundation Ireland under Grant number
18/CRT/6222

Zachary Dair 11 0000-0003-0203-0396 Samantha Dockray 22 0000-0002-0804-8362
Ruairi O’Reilly 11 0000-0001-7990-3461

###### Abstract

Psychophysiology investigates the causal relationship of physiological changes
resulting from psychological states. There are significant challenges with
machine learning-based momentary assessments of physiology due to varying data
collection methods, physiological differences, data availability and the
requirement for expertly annotated data. Advances in wearable technology have
significantly increased the scale, sensitivity and accuracy of devices for
recording physiological signals, enabling large-scale unobtrusive
physiological data gathering. This work contributes an empirical evaluation of
signal variances acquired from wearables and their associated impact on the
classification of affective states by (i) assessing differences occurring in
features representative of affective states extracted from electrocardiograms
and photoplethysmography, (ii) investigating the disparity in feature
importance between signals to determine signal-specific features, and (iii)
investigating the disparity in feature importance between affective states to
determine affect-specific features. Results demonstrate that the degree of
feature variance between ECG and PPG in a dataset is reflected in the
classification performance of that dataset. Additionally, beats-per-minute,
inter-beat-interval and breathing rate are identified as common best-
performing features across both signals. Finally feature variance per-
affective state identifies hard-to-distinguish affective states requiring one-
versus-rest or additional features to enable accurate classification.

###### Keywords:

Machine Learning Classification Psychophysiology Electrocardiogram
Photoplethysmography Affective States

##  1 Introduction

A significant goal of Affective Computing is to improve human-to-computer
interaction by providing a system with a level of emotional intelligence that
aids natural communications and is capable of including emotional components
[27]. This has commonly been approached by deriving emotional states from
speech, facial expressions, gestures and body posture analysis. However,
utilising physiological signals to communicate psychological information is a
recent exploration in the domain, likely due to the increased accessibility of
signals from wearables.

A physiological signal represents an individual’s biological processes derived
from core aspects of human biology. These signals can enable diagnostics, for
instance, analysing heart rate (HR) to detect arrhythmia [29]. Psychological
analysis can also be enabled as mental states originating from unconscious
effort typically present a noticeable physiological change in the relevant
human system [16]. The combined analysis enables a richer understanding of
individuals in terms of their mental and physical health [8].

Psychological states are complex processes comprised of several components,
including feelings, cognitive reactions, behaviour and thoughts[1]. Mapping
psychological states to individual experience provides valuable information
regarding well-being, health (physical and mental), social contexts,
experiences and emotions [7].

Electrocardiograms (ECG) are physiological signals that measure the electrical
activity of the heart. Typically recorded in a clinical setting using multiple
electrodes attached to the individual. Photoplethysmography (PPG) is a
physiological signal used to measure heart activity through variations in the
blood volume of the skin, using a light-emitting-diode and photodetector.
Wearable devices predominately utilise PPG to monitor heart activity. However,
recently advanced wearables have included ECG capabilities for a limited
number of commercial off-the-shelf (COTS) devices.

Data variances occur when recording ECG and PPG due to differing sensor
placement and signal granularity [22, 8]. A lower sampling frequency is
commonly used in PPG compared to ECG to reduce battery consumption in COTS
devices. Such variances are under-recognised in the field of psychophysiology.

This work investigates the impact of signal variances occurring in ECG and PPG
signals acquired from wearable devices for classifying affective states by
addressing the following research aims: (i) To assess differences in features
representative of affective states on a per signal basis, (ii) To investigate
the disparity in precedence ordering of feature importance per signal, and
(iii) To investigate the disparity precedence ordering of feature importance
per affective state.

These aims inform the development of machine learning (ML) pipelines for
classifying affective states. Utilising feature variance per signal to
identify abnormal signal activity or similar affective states which are
causing reduced classification accuracy. In conjunction, feature importance is
utilised to provide insights into feature selection, aiding performance in
tailored signal- or affect-specific approaches.

##  2 Related Work

###  2.1 Heart-Related Physiological Signals

The prevalence of heart-related data in wearable devices stems from a desire
to monitor health through arrhythmia detection and HR as a measure of fitness
[29]. As the heart is controlled involuntarily through the autonomic nervous
system (ANS), it facilitates identifying relationships between involuntary
physiological changes in heart activity and psychological states such as
emotions or behaviour. Multiple psychophysiological theories aim to explain
this relationship, such as Polyvagal Theory [30], which proposes that the ANS
provides the neurophysiological substrates for adaptive behavioural strategies
[28].

Heart activity is complex to capture. In medicine, the gold standard utilises
a 12-lead ECG, resulting in comprehensive data recorded from multiple
electrodes on the human body. However, in ambulatory research and daily life,
this method is not feasible. Typically research-grade (RG) equipment uses
several electrodes, commonly 3-lead ECG, and occasionally includes PPG as an
additional measure. COTS devices tend to rely solely on PPG to monitor heart
activity. However, with recent advances, top-of-the-range smart-watches (Apple
Watch 4-9, Galaxy Active 2, Fitbit) include a 1-lead ECG, which is promising
for portable ECG analysis [25].

Additional physiological signals such as electrodermal activity (EDA),
respiration, skin temperature, electromyogram (EMG), and electrooculogram
(EOG) have demonstrated potential for affective state detection [8, 20];
however, due to additional sensor requirements they are excluded from this
work.

Numerous studies of affective states conduct custom data collection, providing
precise control over the psychological domain explored. Varied stimuli have
aided the elicitation of psychological states, for example, images, movie
clips, music, and dedicated tasks to elucidate stress, such as the Trier
Social Stress Test [33, 2]. As denoted in Table 1, several open-access or on-
request datasets containing ECG and PPG are available. The distinct lack of
emotionally labelled ECG signals from COTS devices is likely due to the recent
inclusion of ECG monitoring capabilities [25].

Table 1: Datasets containing affectively labelled ECG, PPG or both Dataset | ECG | PPG | Participants | Label
---|---|---|---|---
CASE[34] | ✓(1000Hz) | ✓(1000Hz) | 30 | Aro./Val.
WESAD[33] | ✓(700Hz) | ✓(64Hz) | 15 | B, S, A, M
DREAMER[15] | ✓(256Hz) | x | 23 | Aro./Val./Dom.
SWELL[18] | ✓(2048Hz) | x | 25 | S, Aro./Val./Dom.
DEAP[17] | x | ✓(256Hz) | 32 | Aro./Val.
B: Baseline, S: Stress, A: Amusement, M: Meditation
Aro: Arousal, Val: Valence, Dom: Dominance
Table 2: Handcrafted, automated and statistical features utilised for
affective state classification. Note the divergence between features used by
PPG and ECG.

| ECG | PPG | Combined
---|---|---|---
| [31] | [1] | [13] | [5] | [26] | [12] | [6] | [32] | [11] | [24] | [10] | [14] | [21] | [19] | [32] | [4] | [9]
HRV Time: |  |  |
HR | x | x | x | ✓ | x | ✓ | x | ✓ | x | x | x | ✓ | x | x | ✓ | x | x
IBI | x | x | x | x | x | x | x | ✓ | x | x | x | x | x | x | ✓ | ✓ | x
RR | x | x | x | ✓ | ✓ | ✓ | ✓ | ✓ | x | x | x | x | x | x | ✓ | x | ✓
SD | x | x | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | x | x | x | x | x | x | ✓ | x | ✓
P-QRS-T | x | x | ✓ | x | x | x | ✓ | ✓ | x | x | x | x | x | x | ✓ | x | x
HRV Frequency: |  |  |
Low Freq. | x | x | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | x | x | x | ✓ | x | x | ✓ | x | ✓
High Freq. | x | x | x | ✓ | ✓ | ✓ | ✓ | ✓ | x | x | x | ✓ | x | x | ✓ | x | ✓
Freq. Ratios | x | x | x | ✓ | ✓ | ✓ | ✓ | ✓ | x | x | x | ✓ | x | x | ✓ | ✓ | ✓
Non-Linear: |  |  |
Poin. | x | x | x | x | ✓ | ✓ | ✓ | x | x | x | ✓ | x | x | x | x | x | x
Entr. | x | x | x | x | ✓ | ✓ | x | ✓ | ✓ | x | x | x | x | x | ✓ | ✓ | ✓
Deep Learning: |  |  |
Model | ✓ | x | x | x | x | x | x | x | x | x | x | x | ✓ | ✓ | x | x | x
AE | x | x | x | x | x | x | x | x | x | ✓ | x | x | x | x | x | x | x
Additional Features: |  |  |
BR | x | x | x | x | x | ✓ | x | x | x | x | x | x | x | x | x | x | ✓
Sig. Amp. | x | x | x | x | x | x | x | ✓ | x | x | x | x | x | x | ✓ | x | ✓
EMD | x | ✓ | x | x | x | x | ✓ | x | x | x | x | x | x | x | x | ✓ | x
Classifier | CNN | LDA | KNN | LDA | QDC | SVMb | Ens. | SVM | PNN | SVM | SVM | DNN | CNN | CNN | SVM | SVMb | FNN
No. Classes | 4 | 2 | 2 | 3 | 4 | 2 | 4 | 2 | 14a | 4a | 2 | 5 | 2 | 2 | 2 | 2 | 4
Accuracy | 95% | 89% | 85% | 85% | 84% | 82% | 80% | 69% | 100% | 99% | 96% | 91% | 83% | 76% | 65% | 93% | 70%
Datasets | [33] | Priv. | Priv. | Priv. | Priv. | Priv. | Priv. | Priv. | [17] | [33] | [17] | Priv. | [33] | [17] | Priv. | Priv. | Priv.
HRV: Heart Rate Variability, IBI: Inter-beat Interval, RR: R-R Intervals, SD:
Successive Differences,
Freq: Frequency, Poin: Poincare, Entr: Entropy, AE: Auto-Encoder, BR:
Breathing Rate, Ens: Ensemble,
Sig. Amp: Signal Amplitude, EMD: Empirical Mode Decomposition, LDA: Linear
Discriminant Analysis,
(C/D/F)NN: Convolutional/Deep/Feed-Forward Neural Network, QDC: Quadratic
Discriminant Classifier,
SVM: Support Vector Machine, KNN: K-Nearest Neighbours, a One Vs Rest, b
Least-Squares SVM

###  2.2 Affective ECG Analysis

ECG signals contain noise introduced by motion artefacts, biological
differences and sensor de-attachment. Signal processing techniques such as
Butterworth Bandpass, Notch filters and Empirical Mode Decomposition (EMD) are
utilised to reduce the signal noise levels [1]. Subsequently, features
suitable for affective state classification can be extracted from the pre-
processed signals.

An overview of features derivable from ECG and PPG is denoted in Table 2,
grouped by extraction method. Performant ECG-based approaches typically
utilise handcrafted features, particularly time-based HRV features, such as
R-R intervals (RR) which are the intervals between heartbeats, successive
differences (SD) and frequency-based features, such as relative, peak and
absolute power of various frequency bands. Automated feature extraction is
less frequently adopted, with only three of the reviewed approaches utilising
deep learning or signal-processing feature extraction methods.

Recent approaches have favoured deep learning methodologies [31], achieving
significant accuracies on multi-class classifications. However, older studies
focusing on linear and quadratic discriminant analysis (LDA, QDA) [1, 26, 5]
and support vector machines (SVM) [12] remain highly relevant, achieving high
accuracy for their respective classifications. Combinations of ML classifiers
forming ensembles have demonstrated potential for binary classifications in
emotion detection [6]. In comparison to other studies, [31] achieved the
highest accuracy for multiple emotion detection from ECG data utilising a CNN
and reported setting the new state of the art for ECG emotion detection.
Despite the high performance of deep learning approaches in the literature,
this work focuses on classifiers using handcrafted features.

###  2.3 Affective PPG Analysis

PPG analysis provided by COTS devices has typically focused on tracking
medical conditions, physical activity, and stress. The detrimental effects of
stress on human health are a significant motivator for physiological analysis
and preventative healthcare research [3]. However, instances of PPG have
demonstrated similar noise levels to ECG, with the addition of skin tone and
environmental light effects impacting signal quality, requiring signal
cleaning techniques.

There is no consensus on the most frequently used features from the reviewed
PPG-based approaches, see Table 2. The most performant approach [11] leverages
handcrafted non-linear entropy features, followed by [24] using an autoencoder
method for automatic feature extraction. Importantly, both handcrafted and
automatically extracted features aid in achieving a high classification
accuracy above 90% [11, 24, 10, 14].

These affective state classifications are conducted by variations of neural
networks [11, 14, 21, 19] and SVMs [24, 10, 32], which demonstrates great
potential for both binary and multi-class affective state detection using PPG
solely. Notably, these approaches leverage extensive signal processing to
reduce signal noise and contribute to the high performances achieved.

##  3 Methodology

The proposed methodology provides an approach for investigating ECG and PPG
variances and the subsequent impact on affective state classification. The
baseline performance of affective state classification is achieved using
multiple ML classifiers per signal. The inter-signal performance variances are
investigated by analysing the disparity in features between temporally aligned
ECG and PPG, where the degree of feature variance is an indicator of signal
quality. Inter-affective state feature variance is analysed using statistical
measures to provide insights into the distribution and similarity of affective
states. Feature importance is employed to identify commonalities among the
best-performing features across both signals and evaluate each feature’s
utility for affect-specific approaches. Finally, a one-versus-rest (OVR)
classification is adopted to improve performance when classifying similar
affective states.

###  3.1 Datasets

For the purposes of this work, the focus was narrowed to RG physiological
signals due to a lack of publicly available data for COTS devices. “The
Dataset of Continuous Affect Annotations and Physical Signals for Emotion
Analysis” (CASE)[34] and “The Wearable Stress and Affect Detection Dataset”
(WESAD)[33], see Table 1, were utilised in this work. The datasets were
selected due to their inclusion of temporally aligned ECG and PPG with
psychological annotations. Additionally, these signals were recorded using RG
devices in a laboratory environment. CASE incorporates Arousal and Valence
annotations, achieved by collecting joystick movement resulting from
emotionally stimulating video clips. WESAD focuses on stress detection with
limited affective states: a baseline state elicited from “neutral reading”,
amusement caused by comedic video clips, a Trier Social Stress Test[2] to
provoke stress, and a meditation stage aimed at “de-exciting” the individual
following the amusement and stress stages.

###  3.2 Pre-Processing

ECG and PPG signals recorded per subject within these datasets span the
duration of the experiment resulting in approx 91/40 minutes for WESAD/CASE.
Each signal is pre-processed into 10-second windows to facilitate analysis,
accomplished using a sliding window technique with a 1-second overlap. A
10-second duration was selected due to efficient performance demonstrated in
[31]; additionally, this duration enables low latency as classification occurs
every 10-seconds and contains adequate data for feature computation.

A Butterworth-Bandpass filter is used to reduce signal noise, facilitating the
extraction of selected features while maintaining a degree of “rawness” in the
signal. This filter was adopted as it is frequently adopted in the literature
and more closely aligns with COTS devices and their reduced computational
power.

Once filtered and windowed, the data is aligned with the psychological
annotations. For WESAD, annotations were numeric values sampled at 700Hz. Each
value from 0-4 is associated with the psychological states: Transient,
Baseline, Stress, Amusement and Meditation. Annotations 5-7 and Transient data
are omitted as per the author’s instructions [33]. Certain windows may include
multiple emotive annotations; hence to identify the most pertinent emotion,
the mean of all annotation values per window is calculated and rounded to the
nearest annotation (1-4) using Euclidean distance. Alternative approaches[6]
omit these windows and the neighbouring segments to prevent overlap.

A similar procedure is required for CASE; the raw annotation data is provided
as values on an x and y-axis representing Arousal and Valence[34], these
values are normalised to a range of 0.5 to 9.5, and subsequently converted to
discrete representations, resulting in low (0.5-5) and high (5.01-9.5) Arousal
and Valence for each window.

Both signals provide capabilities to derive a wide array of handcrafted
features useful for identifying affective information. This work utilises a
python toolkit HeartPy [35] to enable extraction of HRV features from each
window of data, summarised in accompanying table of Figure 4.

###  3.3 HRV Feature Variance

The feature variance approach proposed is to statistically evaluate any
disparity occurring in derived features from ECG and PPG under multiple
conditions. Inter-signal variance is evaluated by computing the absolute
difference between an ECG-derived feature and its PPG counterpart from
temporally aligned signals. This is assessed using the same window of heart
activity and provides a granular analysis to aid in identifying noisy, erratic
or abnormal signal activity, causing unreliable computations of features. This
variance is depicted by a significant absolute difference of a feature between
the two signals.

Analysing the inter-affective state variance in features enables the
identification of the degree of change between states, as investigated in [10,
12]. The proposed methodology computes the minimum, maximum, mean and standard
deviation of each feature value per affective state. Additionally outliers are
identified, these are observations found in the upper and lower quartiles.
This method identifies states which are complex to distinguish due to a
similar feature distribution, such as meditation and relaxation. An OVR
approach is adopted to convert a multi-class problem into multiple binary
classifications. Using OVR, a classifier aims to identify an affective state
individually from the remaining states, which increases the degree of
distinction between classes.

###  3.4 HRV Feature Importance

This work adopts a game theory approach for feature importance known as
“Shapley Additive exPlanations“[23]. This method computes SHAP values
representing the degree of change on the classifier output caused by each
individual feature, the magnitude of change and number of samples affected
indicate the impact factor of a given feature.

Feature importance has enabled the identification of signal-specific features
in [32]. However, their approach utilised different features for ECG and PPG,
as such, an intra-signal comparison could not be conducted, which would
provide insights into the commonality of features between ECG and PPG,
motivating the intra-signal feature importance analysis provided in this work.

Feature importance can also provide insights into the variance of features per
affective state, valuable for the creation of tailored emotion-specific
approaches. In [9], a neural network is used for classification, and the most
important features were identified from the first layer’s weights. These
features were then evaluated to identify a statistical difference between
affective states.

###  3.5 ML Based Classification of Affective State

A range of classifiers was selected to provide a holistic view of the
classification performance using different architectures and the suitability
of ECG and PPG for automated affective state detection.

Each classifier conducts a per-signal classification on each dataset, where
20% of the data acts as a hold-out test set, which is unseen data used to
evaluate the final classifier. To ensure generalisability, five-fold cross-
validation is utilised, transforming the remaining 80% of data into “folds”,
enabling a per-fold classification. Subsequently, comparing the per-fold and
average performance across the five folds enables the identification of the
most robust and performant classifier. Finally, the most performant classifier
is trained on the entire training set and evaluated against the hold-out set
to assess expected performance in real-life classifications.

##  4 Results and Discussion

Figure 1: Mean cross-validation accuracy classifying affective states for
model selection, and the performance on the holdout test set from the best
performing classifier. The performance variance in ExtraTrees classifier
between ECG and PPG reflects the degree of feature variance identified per
dataset.

###  4.1 HRV Feature Variance

Figure 2: Absolute difference between ECG and PPG features: BPM, IBI and BR,
for CASE and WESAD. High variance demonstrates unreliable feature computation
in one of the signals due to signal noise or sensor differences.

(a) Outliers Included

(b) Outliers Excluded

Figure 3: Inter-signal and inter-affective state variance for BPM in WESAD,
including and excluding outliers. Note in (a) the presence of outliers with a
BPM of over 300 occurring in ECG indicating abnormal signal activity.
Additionally, in (a, b), a visible overlap in neutral, amusement, and
meditation occur, demonstrating the degree of similarity in these states.

The wearables’ sample rate disparity (See Table 1) is evident in the inter-
signal feature variance results depicted in Figure 2. The reduced sample rates
in WESAD result in slightly decreased granularity of ECG data and
significantly in PPG data compared to CASE. A higher fluctuation in feature
variance occurs in WESAD in terms of magnitude and frequency, stemming from
the high sample rate disparity.

In CASE, beats per-minute (BPM) and inter-beat interval (IBI) contain a small
variance with substantial spikes relative to the average. These variances
occur in isolated data segments and are likely caused by electrode
disconnection, movement, or subject-specific factors, visible in Figure 2 at
approximately window numbers 250, 550 and 900. Such occurrences may benefit
from additional signal processing to reduce noise and improve feature
computation accuracy.

Interestingly, the breathing rate (BR) feature exhibits a high deviation
between signals in both datasets. This deviation indicates that at least one
of the signals is unreliably computing BR, likely due to the wrist and finger
placement of the PPG sensors.

A low degree of inter-affective state feature variance was identified between
WESAD baseline, amusement, and meditation states for all features indicating
these states are difficult to distinguish as depicted in Figure 3.
Statistically similar features negatively impact automated classification, as
the classifier struggles to differentiate between the classes. This impact is
demonstrated by the reduced performance in multi-class classifications
(58%-69%) as compared to the OVR performance depicted by the ROC curves (ROC
Area: 0.70-0.95) in Figure 5. This performance increase validates the utility
of OVR classifications when classifying affective states that are difficult to
differentiate due to statistical similarities.

###  4.2 HRV Feature Importance

Analysing the SHAP values per feature indicates that BPM, IBI, and BR have the
most significant impact on classification for both signals, as demonstrated in
Figure 4. The remaining features exhibit inconsistent influence between the
signals. Most notably, standard deviation 1 divided by standard deviation 2
(SD1/SD2) and room-mean-square of successive differences (RMSSD) exhibit
higher impact in PPG as opposed to ECG. This demonstrates the need for
assessing feature importance on a per-signal basis to identify which features
are most informative for use in tailored signal-specific classification
approaches.

Certain features demonstrate varying impacts across affective states,
indicating the presence of affect-specific features. For example, BPM and IBI
exhibit high impacts on the class “stress”, indicating their suitability for
stress detection approaches. Assessing feature importance per-affective state
provides an informative analysis of feature utility for affect-specific
approaches.

The high feature importance of BPM for “stress” is due to statistical
distinction to the other affective states in the inter-affective state feature
variance, as depicted in Figure 3. This demonstrates the benefit of assessing
inter-affective state feature variance and feature importance to gain insights
to aid the creation of affect-specific approaches.

Feature | Abbrv.
---|---
Beats Per Minute | BPM
Interbeat Interval | IBI
Std dev. of RR Intervals | SDNN
Std dev. of successive diff. | SDSD
RMSE of successive diff. | RMSSD
Proportion of diff. $<$ 20ms | pNN20
Proportion of diff. $<$ 50ms | pNN50
Median absolute dev. of RR intervals | MAD
Estimated breathing rate | BR
Poincare analysis | SD1, SD2,
| S, SD1/SD2
Dev: Deviation, Diff: Difference |

Figure 4: Shapley Additive exPlanations (SHAP) Feature Importance from
ExtraTrees classifying WESAD signals.

###  4.3 Automated Affective State Classification Variance

Finally, the selected classifier is trained on the initial 80% of data and
classifies the hold-out set to assess expected performance in real-life
classifications. The ExtraTrees classifier (ET) was selected as the most
performant classifier from the model selection, where it was trained on 80% of
the training data and evaluated on the remaining 20%. Notably, ET exhibits an
increased performance when evaluated on the hold-out set as it was trained on
all available training data. The full model comparison and ET hold-out
performance is depicted in Figure 1. Interestingly, the classifier performance
variance between ECG and PPG is similar to the degree of the inter-signal
feature variance identified per dataset.

In contrast with the state-of-the-art [31, 11], the performance achieved is
lower for ECG and PPG; however, this work focuses on the analysis and
understanding of variances between the signals for affective analysis rather
than achieving high classification accuracy. Analysing the ROC curves from ET
demonstrates the true and false positive rates per signal for each affective
state, see Figure 5. On average, ECG demonstrates increased capabilities for
affective classification by achieving a higher ROC area than PPG, varying with
a range of 0.02-0.11. The increased performance via OVR demonstrates the
benefit of identifying and overcoming the effects of similar affective states
to achieve greater classification performance.

(a) WESAD

(b) CASE

Figure 5: ROC Curves from ExtraTrees representing the OVR classification
variance between ECG and PPG

##  5 Conclusions

The inter-signal classification performance disparity mirrors the degree of
feature variance between signals from both datasets. Specifically, WESAD
exhibited a high feature variance, which explains the higher disparity in
classification accuracy and ROC area per signal. Conversely, a lower inter-
signal feature variance and a lower disparity in the performance measures
occurred for CASE. This demonstrates the utility of inter-signal feature
variance in identifying inconsistent computations of features stemming from
sensor differences or abnormal signal activity, which negatively impact
classification performance. These occurrences are likely to be more frequent
in the ambulatory analysis due to motion artefacts and uncontrolled usage of
wearables.

Furthermore, inter-affective state feature variance enables the identification
of affective states that contain a similar distribution of features, which
causes classification confusion. To counter this, the similar states are
aggregated into an OVR classification problem, leading to increased
performance, demonstrated by the ROC area per affective state.

Feature importance identifies BPM, IBI, and BR as the most impactful features
for affective classification across ECG and PPG. Notably, the remaining
features exhibit inconsistent impacts, specifically SD1/SD2 and RMSSD, which
demonstrate a greater impact in PPG, warranting the exploration of signal-
specific features. Analysing statistical measures to understand the inter-
affective state feature variance indicates that certain features provide a
greater degree of affect-specific information beneficial for tailored
applications.

This work contributes an empirical analysis of data variances in ECG and PPG
acquired using wearables and the impact on affective state classification.
Therefore, enabling practitioners to make informed decisions when creating ML
pipelines for affective state classification. The code-base will be made open
access on Github (https://github.com/ZacDair/Emo-Phys-Eval), enabling
automated feature variance analysis from each of these perspectives in a
combined manner, regardless of data acquisition methods. While this approach
analyses handcrafted features, it can also be utilised with automatically
extracted features.

Future work will expand the analysis by utilising additional datasets to
provide greater insights into the variances stemming from data collection
devices, affective states, and population differences. In addition, an
extended analysis will be conducted using additional features and methods to
further inform the development of ML pipelines for affective state detection.

## References

  * [1] Agrafioti, F., et al.: Ecg pattern analysis for emotion detection. IEEE Trans. Affective Comput. 3(1), 102–115 (2012). https://doi.org/10.1109/T-AFFC.2011.28
  * [2] Birkett, M.A.: The trier social stress test protocol for inducing psychological stress. J. Vis. Exp. (Oct 2011). https://doi.org/10.3791/3238
  * [3] Can, Y.S., Chalabianloo, N., Ekiz, D., Ersoy, C.: Continuous stress detection using wearable sensors in real life: Algorithmic programming contest case study. Sensors 19(8), 1849 (2019). https://doi.org/10.3390/s19081849
  * [4] Cheema, A., Singh, M.: Psychological stress detection using phonocardiography signal: An empirical mode decomposition approach. Biomed. Signal Process. Control 49, 493–505 (2019). https://doi.org/10.1016/j.bspc.2018.12.028
  * [5] Cinaz, B., Arnrich, B., La Marca, R., Tröster, G.: Monitoring of mental workload levels during an everyday life office-work scenario. Personal and Ubiquitous Computing 17 (2013). https://doi.org/10.1007/s00779-011-0466-1
  * [6] Dissanayake, T., Rajapaksha, Y., Ragel, R., Nawinne, I.: An ensemble learning approach for electrocardiogram sensor based human emotion recognition. Sensors 19(20) (2019). https://doi.org/10.3390/s19204495
  * [7] Dockray, S., O’Neill, S., Jump, O.: Measuring the psychobiological correlates of daily experience in adolescents. Journal of Research on Adolescence 29(3), 595–612 (2019). https://doi.org/10.1111/jora.12473
  * [8] Dzedzickis, A., Kaklauskas, A., Bucinskas, V.: Human emotion recognition: Review of sensors and methods. Sensors 20(3) (2020). https://doi.org/10.3390/s20030592
  * [9] Filippini, C., Di Crosta, A., Palumbo, R., Perpetuini, D., Cardone, D., Ceccato, I., Di Domenico, A., Merla, A.: Automated affective computing based on bio-signals analysis and deep learning approach. Sensors 22(5) (2022). https://doi.org/10.3390/s22051789
  * [10] Goshvarpour, A., Goshvarpour, A.: Poincaré’s section analysis for ppg-based automatic emotion recognition. Chaos, Solitons & Fractals 114, 400–407 (2018). https://doi.org/10.1016/j.chaos.2018.07.035
  * [11] Goshvarpour, A., Goshvarpour, A.: Evaluation of novel entropy-based complex wavelet sub-bands measures of ppg in an emotion recognition system. Journal of Medical and Biological Engineering 40(3), 451–461 (Jun 2020). https://doi.org/10.1007/s40846-020-00526-7
  * [12] Hsu, Y.L., Wang, J.S., Chiang, W.C., Hung, C.H.: Automatic ECG-based emotion recognition in music listening. IEEE Trans. Affective Comput. 11(1), 85–99 (2020). https://doi.org/10.1109/TAFFC.2017.2781732
  * [13] Jing, C., Liu, G., Hao, M.: The research on emotion recognition from ecg signal. In: 2009 International Conference on Information Technology and Computer Science. vol. 1, pp. 497–500 (2009). https://doi.org/10.1109/ITCS.2009.108
  * [14] Kalra, P., Sharma, V.: Mental stress assessment using ppg signal a deep neural network approach. IETE Journal of Research 0(0), 1–7 (2020). https://doi.org/10.1080/03772063.2020.1844068
  * [15] Katsigiannis, S., Ramzan, N.: Dreamer: A database for emotion recognition through eeg and ecg signals from wireless low-cost off-the-shelf devices. IEEE Journal of Bio. and Health Info. (2018). https://doi.org/10.1109/JBHI.2017.2688239
  * [16] Kim, J., André, E.: Emotion recognition based on physiological changes in music listening. IEEE PAMI 30(12), 2067–2083 (2008). https://doi.org/10.1109/TPAMI.2008.26
  * [17] Koelstra, S., Muhl, C., Soleymani, M., Lee, J.S., Yazdani, A., Ebrahimi, T., Pun, T., Nijholt, A., Patras, I.: Deap: A database for emotion analysis ;using physiological signals. IEEE Trans. Affective Comput. 3(1), 18–31 (2012). https://doi.org/10.1109/T-AFFC.2011.15
  * [18] Koldijk, S., Sappelli, M., Verberne, S., Neerincx, M.A., Kraaij, W.: The SWELL knowledge work dataset for stress and user modeling research. In: ICMI. pp. 291–298. ACM (2014). https://doi.org/10.1145/2663204.2663257
  * [19] Lee, M.S., Lee, Y.K., Pae, D.S., Lim, M.T., Kim, D.W., Kang, T.K.: Fast emotion recognition based on single pulse ppg signal with convolutional neural network. Applied Sciences 9(16) (2019). https://doi.org/10.3390/app9163355
  * [20] Lin, S., et al.: A review of emotion recognition using physiological signals. Sensors 18(7), 2074 (Jun 2018). https://doi.org/10.3390/s18072074
  * [21] Lisowska, A., et al.: Catching patient’s attention at the right time to help them undergo behavioural change. In: Artificial Intelligence in Medicine, pp. 72–82. Spr. Int. Pub. (2021), https://link.springer.com/10.1007/978-3-030-77211-6_8
  * [22] Mahdiani, S., et al.: Is 50 hz high enough ecg sampling frequency for accurate hrv analysis? In: EMBC. pp. 5948–5951 (2015). https://doi.org/10.1109/EMBC.2015.7319746
  * [23] Molnar, C.: Interpretable Machine Learning. Lulu.com (2022), christophm.github.io/interpretable-ml-book/
  * [24] Mukherjee, N., et al.: Real-time mental stress detection technique using neural networks towards a wearable health monitor. Meas. Sci. Technol. 33(4), 044003 (2022). https://doi.org/10.1088/1361-6501/ac3aae
  * [25] Nabeel, S., et al.: A comparison of manual electrocardiographic interval and waveform analysis in lead 1 of 12-lead ecg and apple watch ecg: A validation study. Cardiovascular Digital Health Journal (2020). https://doi.org/10.1016/j.cvdhj.2020.07.002
  * [26] Nardelli, M., Valenza, G., Greco, A., Lanata, A., Scilingo, P.: Recognizing emotions induced by affective sounds through heart rate variability. IEEE Trans. Affective Comput. 6(4), 385–394 (2015). https://doi.org/10.1109/TAFFC.2015.2432810
  * [27] Picard, R.W.: Affective computing: challenges. Int. J. Hum. Comput. Stud. 59(1), 55–64 (2003). https://doi.org/10.1016/S1071-5819(03)00052-1
  * [28] Porges, S.W.: The polyvagal theory: new insights into adaptive reactions of the autonomic nervous system. Clevel. Clin. J. Med. 76 Suppl 2(Suppl 2), S86–S90 (Apr 2009). https://doi.org/10.3949/ccjm.76.s2.17
  * [29] José da S. Luz, E., et al.: Ecg-based heartbeat classification for arrhythmia detection: A survey. CMPB 127, 144–164 (2016). https://doi.org/10.1016/j.cmpb.2015.12.008
  * [30] S W, P., et al.: Vagal tone and the physiological regulation of emotion. Monogr. Soc. Res. Child Dev. 59(2-3), 167–186 (1994)
  * [31] Sarkar, P., Etemad, A.: Self-supervised ECG representation learning for emotion recognition. IEEE Trans. Affective Comput. (2021). https://doi.org/10.1109/TAFFC.2020.3014842
  * [32] Sayed Ismail, S.N.M., Ab. Aziz, N.A., Ibrahim, S.Z.: A comparison of emotion recognition system using electrocardiogram (ecg) and photoplethysmogram (ppg). J. King Saud Univ. - Comput. Inf. Sci. 34(6, Part B), 3539–3558 (2022). https://doi.org/10.1016/j.jksuci.2022.04.012
  * [33] Schmidt, P.a.a.: Introducing wesad, a multimodal dataset for wearable stress and affect detection. In: ICMI 20. ACM (2018). https://doi.org/10.1145/3242969.3242985
  * [34] Sharma, K., et al.: A dataset of continuous affect annotations and physiological signals for emotion analysis (2018). https://doi.org/10.48550/ARXIV.1812.02782
  * [35] van Gent, P., et al.: Heartpy: A novel heart rate algorithm for the analysis of noisy signals. Transp. Res. F: Traffic Psychol. Behav. 66, 368–378 (2019). https://doi.org/10.1016/j.trf.2019.09.015


# Discovery of two Einstein crosses from massive post–blue nugget galaxies at
$z>1$ in KiDS111Based on observations with OmegaCam@VST and MUSE@VLT (Prog.
ID: 0105.A-0253).

N.R. Napolitano School of Physics and Astronomy, Sun Yat-sen University Zhuhai
Campus, Daxue Road 2, 519082—Tangjia, Zhuhai, Guangdong, China
INAF—Osservatorio Astronomico di Capodimonte, Salita Moiariello 16,
I-80131—Napoli, Italy  R. Li  School of Physics and Astronomy, Sun Yat-sen
University Zhuhai Campus, Daxue Road 2, 519082—Tangjia, Zhuhai, Guangdong,
China  C. Spiniello  Department of Physics, University of Oxford, Denys
Wilkinson Building, Keble Road, Oxford OX1 3RH, UK  INAF—Osservatorio
Astronomico di Capodimonte, Salita Moiariello 16, I-80131—Napoli, Italy  C.
Tortora INAF—Osservatorio Astronomico di Capodimonte, Salita Moiariello 16,
I-80131—Napoli, Italy  INAF—Osservatorio Astronomico di Arcetri, L.go E. Fermi
5, I-50125—Firenze, Italy  A. Sergeyev  Institute of Astronomy, V. N. Karazin
Kharkiv National University, 35 Sumska Str., Kharkiv, Ukraine  Institute of
Radio Astronomy of the National Academy of Sciences of Ukraine  G. D’Ago
Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña
Mackenna 4860, 7820436 Macul, Santiago, Chile  X. Guo School of Astronomy and
Space Science, Nanjing University, Nanjing, Jiangsu 210093, China  Key
Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry
of Education, Nanjing 210093, China  L. Xie  School of Physics and Astronomy,
Sun Yat-sen University Zhuhai Campus, Daxue Road 2, 519082—Tangjia, Zhuhai,
Guangdong, China  M. Radovich  INAF - Osservatorio Astronomico di Padova,
vicolo dell’Osservatorio 5, I-35122 Padova, Italy  N. Roy  School of Physics
and Astronomy, Sun Yat-sen University Zhuhai Campus, Daxue Road 2,
519082—Tangjia, Zhuhai, Guangdong, China  L. V. E. Koopmans  Kapteyn
Astronomical Institute, University of Groningen, P.O.Box 800, 9700AV
Groningen, the Netherlands  K. Kuijken  Leiden Observatory, Leiden University,
P.O.Box 9513, 2300RA Leiden, The Netherlands  M. Bilicki Center for
Theoretical Physics, Polish Academy of Sciences, al. Lotników 32/46, 02-668,
Warsaw, Poland  T. Erben  Argelander-Institut für Astronomie, Auf dem Hügel
71, 53121 Bonn / Germany  F. Getman  INAF—Osservatorio Astronomico di
Capodimonte, Salita Moiariello 16, I-80131—Napoli, Italy  C. Heymans
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford
Hill, Edinburgh, EH9 3HJ, UK  Ruhr-University Bochum, Astronomical Institute,
German Centre for Cosmological Lensing, Universitätsstr. 150, 44801 Bochum,
Germany  H. Hildebrandt  Ruhr-University Bochum, Astronomical Institute,
German Centre for Cosmological Lensing, Universitätsstr. 150, 44801 Bochum,
Germany  C. Moya  Instituto de Astrofísica, Pontificia Universidad Católica de
Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile  H.Y. Shan
Shanghai Astronomical Observatory (SHAO), Nandan Road 80, Shanghai 200030,
China  University of Chinese Academy of Sciences, Beijing 100049, China  G.
Vernardos  Institute of Astrophysics, Foundation for Research and
Technology$-$ Hellas (FORTH), GR-70013, Heraklion, Greece  A.H. Wright  Ruhr-
University Bochum, Astronomical Institute, German Centre for Cosmological
Lensing, Universitätsstr. 150, 44801 Bochum, Germany

(Accepted ApJL)

###### Abstract

We report the discovery of two Einstein Crosses (ECs) in the footprint of the
Kilo-Degree Survey (KiDS): KIDS J232940-340922 and KIDS J122456+005048. Using
integral field spectroscopy from MUSE@VLT, we confirm their gravitational-lens
nature. In both cases, the four spectra of the source clearly show a
prominence of absorption features, hence revealing an evolved stellar
population with little star formation. The lensing model of the two systems,
assuming a singular isothermal ellipsoid (SIE) with external shear, shows
that: 1) the two crosses, located at redshift $z=0.38$ and 0.24, have Einstein
radius $R_{\rm E}=5.2$ kpc and 5.4 kpc, respectively; 2) their projected dark
matter fractions inside the half effective radius are 0.60 and 0.56 (Chabrier
IMF); 3) the sources are ultra-compact galaxies, $R_{\rm e}\sim 0.9$ kpc (at
redshift $z_{\rm s}=1.59$) and $R_{\rm e}\sim 0.5$ kpc ($z_{\rm s}=1.10$),
respectively. These results are unaffected by the underlying mass density
assumption. Due to size, blue color and absorption-dominated spectra,
corroborated by low specific star-formation rates derived from optical-NIR
spectral energy distribution fitting, we argue that the two lensed sources in
these ECs are blue nuggets migrating toward their quenching phase.

editorials, notices — miscellaneous — catalogs — surveys

##  1 Introduction

Strong gravitational lensing is a powerful tool to measure the distribution of
Dark Matter (DM) in galaxies and study the properties of high-redshift
sources. According to General Relativity, matter in the universe acts as a
cosmic telescope deforming and magnifying the light of objects which would be
hardly observable otherwise. Depending on the size of the source and the
alignment of the foreground galaxy (lens or deflector) and the source, strong
lensing events show-up as arcs or rings (when the source is extended, e.g. a
galaxy) or as multiple images (when the source is a compact system, e.g. a
quasar). Deformed images of background galaxies can be used, in combination
with the dynamical and stellar population analysis of the deflector, to
determine the lens total mass density profiles (e.g., Koopmans et al. 2006,
Auger et al. 2010, Nightingale et al. 2019), to separate the dark from the
luminous matter and to constrain the lens stellar initial-mass-function (IMF)
slope (e.g., Treu et al. 2010 , Spiniello et al. 2011). Furthermore, doubly
(‘doublets’) or even quadruply-lensed (‘quads’) QSOs are particularly valuable
for cosmology because they enable measuring the Hubble constant ($H_{0}$, Suyu
et al. 2013) via time-delays among the QSO light curves.

However, doublets and quads can also be produced by other compact sources,
e.g. high–redshift, starforming, ultra–compact galaxies (e.g. Muzzin et al.
2012), which are fairly common at high–$z$. Using lensing forecasts from
typical ground based surveys (Collett, 2015) with a depth of the order of
$r\sim 25$, the number of expected quads from compact galaxies ($R_{\rm e}<1$
kpc) at redshift $z\;\hbox to0.0pt{\lower 2.5pt\hbox{$\sim$}\hss}\raise
1.5pt\hbox{$<$}\;4$, with an alignment sufficient to make a cross-like
geometry (e.g. source misalignment $<$0.1′′), is of the order of half a dozen
every 1000 deg2 (see also §4). Unfortunately, only few of such systems have
been observed so far: besides the Muzzin et al. (2012) system, only two
Einstein Cross (EC) configurations from Ly-$\alpha$ emitters have been
confirmed and fully analysed (Bolton et al. 2006, Bettoni et al. 2019). ECs
are interesting per se as they are the rarest and most spectacular
manifestation of quad systems, showing a distinctive symmetric cross pattern
around the deflector, generated when the source and the lens are almost
perfectly aligned. Generally, these systems have been found to be produced by
distant quasars (e.g. Magain et al. 1988, Ostrovski et al. 2018) or supernovae
(aka Refsdal system, Kelly et al. 2015).

Within the Kilo-Degree Survey (KiDS, de Jong et al. 2015), we have undertaken
a systematic search for strong gravitational lenses, both arcs (Petrillo et
al. 2017, P+17 hereafter, Petrillo et al. 2019a, b) and multiple images
(Spiniello et al., 2018, 2019). In particular, in the process of improving the
overall efficiency of the Convolutional Neural Network (CNN) finders, started
with P+17, we have collected, in an area of $\sim$1000deg2, a series of high
quality quad candidates (Li et al., 2020, Li+20 hereafter), among which we
found three clear EC configurations.

In this paper we report the results of the spectroscopic follow-up of the best
two of them. We argue that these represent a new class of sources of EC
configurations, i.e. high$-z$ post-blue nugget systems, and discuss the
possibility to systematically search for these objects in current and future
ground and space surveys. For all calculations, we assume a $\Lambda$CDM
cosmology with ($\Omega_{M}$, $\Omega_{\Lambda}$, $h$)=(0.3, 0.7, 0.7).

##  2 Confirmation and Lensing model

The first EC, KIDS J232940-340922 (KIDS-EC1, hereafter, see Fig. 1), has been
found in the Southern KiDS patch. The deflector has a AB magnitude of $r\sim
19.8$ and red color, $g-i=1.4$, while the average magnitude of the 4 lensed
images of the source is $r\sim 22.6$ and their average color is $g-i=-0.21$.
The second EC, KIDS J122456+005048 (KIDS-EC2, hereafter, Fig. 1), has been
found in the Northern KiDS patch. The deflector has total magnitude $r\sim
19.7$ and color $g-i=1.9$, while the 4 images have an average magnitude $r\sim
22.0$ and average color $g-i=1.0$, i.e. bluer than the deflector but redder
than KIDS-EC1 source. The two sources have obtained a high CNN probability and
also high human visual score (see Li+20), hence they have been selected for
the spectroscopic follow-up. Table 1 lists coordinates, relative positions of
lens and source images and the optical and near-infrared (NIR) photometry in
the $ugriZYJHK_{\rm s}$ bands, for both ECs. To minimize the relative
contamination and derive homogeneous photometry for all sources, the 9-band
photometry of lens and lensed images are derived by a simultaneous seeing
convolved Sérsic plus 2D Gaussian fit of the objects, respectively, from KiDS-
DR4 (Kuijken et al., 2019) and VIKING (Edge et al., 2013) calibrated images.

###  2.1 MUSE spectroscopy and lensing confirmation

Spectroscopic observations have been collected under ESO Directory
Discretionary Time (program ID: 0105.A-0253, PI Napolitano) with MUSE at VLT,
Cerro Paranal. Run A has been completed in November 2019 for KIDS
J232940-340922 and Run B on February 2020 for KIDS J122456+005048.
Observations have been taken in service mode, in wide-field non adaptive optic
configuration, which allows a full $1^{\prime}\times 1^{\prime}$ field of
view, in the wavelength range $\lambda=[4750,9300]$Å. The MUSE grating
spectral resolution varies from 1750 to 3750, end-to-end, in the same
interval. The total exposure time for both targets is 130 minutes divided in 3
observing blocks (OBs). Every OB is split in 2$\times 1300s$ exposures with a
90deg position angle offset. The final seeing of the combined exposures is
0.93′′ for KiDS-EC1 and 0.84′′ for KiDS-EC2.

Reduced data have been provided by ESO as Internal Data Products, using the
official MUSE pipeline (v2.8). For KiDS-EC1 we have also performed our own
data-reduction to check consistency with the ESO Phase 3 data products, using
the same pipeline. We have found very consistent spectral quality, in terms of
flat fielding, signal-to-noise-ratio (SNR) and sky subtraction. A zoom-in of
the field-of-views (FOVs) of the EC datacubes (i.e. the integrated flux over
all wavelengths) are reported in Fig. 1, and clearly show all sources seen in
the KiDS $gri$ color images (also shown in the same figure). In Fig. 1, we
also present the 1D spectra from the four lensed images, extracted from a
single MUSE pixel, and the one of the corresponding deflector, extracted over
an aperture of 3 pixels for KIDS-EC1 and 4 pixels for KIDS-EC2 respectively,
corresponding to about half of the effective radius, $R_{\rm e}/2$, see §2.2.
We use these apertures because they allow us to minimise the contamination
from the blue lensed sources. For both systems, the lensing nature is
confirmed by the presence of identical spectral features in the four different
images, all consistent with the same redshift, higher than the one of the
lens. For KIDS-EC1, we infer a redshift of $z=1.590\pm 0.001$ from Fe and Mg
absorption lines and faint [FeII] ($\lambda=5650$Å), [CII] ($\lambda=6025$Å)
and [OII] ($\lambda=6400$Å) emissions. For KIDS-EC2, we calculate instead
$z=1.102\pm 0.001$, from the spectrum at $\lambda>7800$Å, including Balmer
absorption lines (H10 at $\lambda\sim 7980$Å, H9 at $\lambda\sim 8060$Å, H8 at
$\lambda\sim 8180$Å, H$\epsilon$ at $\lambda\sim 8350$Å and H$\delta$ at
$\lambda\sim 8620$Å), the Calcium K and H doublet (at $\lambda\sim
8360-8450$Å), although with some sky contamination in H, and a clear [OII]
doublet emission line at $\lambda\sim 7850$Å. In this latter case Fe and Mg
absorption lines are also present but look slightly blue-shifted, possibly due
to some gas outflow from the source galaxy (see e.g. Rubin et al. 2014;
Burchett et al. 2020), which we will investigate in a forthcoming paper.

The two deflectors reveal typical features of early-type galaxies (ETGs), in
particular a strong rest-frame break at 4000Å, faint Balmer lines in
absorption, and Fe, Mgb and NaD lines, all characteristic of an old, metal
rich stellar population. The inferred deflector redshifts are $z_{\rm
l}=0.381\pm 0.001$ for KIDS-EC1, and $z_{\rm l}=0.237\pm 0.001$ for KIDS-EC2.
The SNR and resolution of the spectra allow us to estimate the velocity
dispersion of these two systems using the pPXF software (Cappellari, 2017),
yielding $\sigma_{R_{\rm e}/2}=192\pm 4$ km s-1 for KIDS-EC1 and
$\sigma_{R_{\rm e}/2}=248\pm 2$ km s-1 for KIDS-EC2222We estimate that
systematic errors from template mismatch and masked regions may amount to
$\sim 20$ km s-1. (see also Table 2). The best-fit models are overlaid on the
galaxy spectra in the same Fig. 1.

###  2.2 Lensing model, dynamical masses, dark matter fractions

The two ECs are modelled using the lfit$\\_$gui code (Shu et al., 2016). We
use KiDS $r-$band images, with a pixel scale of $0.2^{\prime\prime}$ and
seeing $0.8^{\prime\prime}$ for KIDS-EC1 and $0.7^{\prime\prime}$ for KIDS-
EC2. The effect of the seeing is taken into account by convolving the lensing
models with a point spread function (PSF) generated by nearby stars (see e.g.,
Roy et al., 2018). The lfit$\\_$gui code simultaneously models the deflector
light, the lensed image positions and their magnification, and the best
position and light distribution of the source galaxy. We assume a single
Sérsic profile (Sérsic, 1963) for the two deflectors, although for KIDS-EC2 we
need to account also for the presence of two nearby galaxies (G1 and G2,
respectively, see Fig. 1), as well as for the two sources. The deflector total
mass distributions are modelled with a singular isothermal ellipsoid (SIE,
Kormann et al., 1994) profile with a projected two-dimensional surface mass
density profile described by:

$\Sigma(x,y)=\frac{1}{2}\Sigma_{\rm c}\,\sqrt{q}\,\theta_{\rm
E}\,(x^{2}+q^{2}y^{2})^{-1/2},$ (1)

where $\theta_{E}$ is the lensing strength, equivalent to the Einstein radius,
$q$ is the minor-to-major axis ratio of the isodensity contours. $\Sigma_{\rm
c}=c^{2}/(4\pi G)(D_{\rm s}/D_{\rm d}D_{\rm ds})$ is the critical density,
where $D_{\rm s}$ and $D_{\rm d}$ are the angular diameter distances from the
observer of the lens and the source, respectively, and $D_{\rm ds}$, the
distance between deflector and source. The assumption of a SIE model is
motivated by evidence pointing toward a logarithmic mass-density slope close
to $-2$ for the total mass density around the Einstein radius (e.g. Koopmans
et al. 2006), however the impact of this assumption is discussed in Appendix,
while a full modeling with more general density profiles will be presented in
future detailed analyses. We also include external shear, $\gamma_{\rm ext}$,
which approximates the influence of the surrounding environment on the lensing
potential. After having initialized the model, with some test runs, the final
best-fit of the two ECs are obtained via $\chi^{2}$ minimization using the
Levenberg–Marquardt algorithm (Moré, 1978). These are shown in Fig. 2, and the
corresponding parameters are reported in Table 2. From the Table, we can draw
some first general results: 1) the two ECs show similar Einstein radii,
$R_{\rm E}$, which are both very close to the lens effective radii, $R_{\rm
e}$ (i.e. $R_{\rm E}/R_{\rm e}\sim 0.9$); 2) the stellar velocity dispersion
measured from the spectrum of KIDS-EC1 is smaller than that inferred by the
lens model, suggesting that the actual slope of the total density profile
might deviate from $-2$ (see e.g. Auger et al. 2010 and the discussion in the
Appendix), while the two values are consistent within 2$\sigma$ for KIDS-EC2;
3) for both lenses, the total mass may be rounder than the starlight
distribution (i.e. $b/a$=0.91 vs. 0.89 for KIDS-EC1 and $b/a$=0.68 vs. 0.59
for KIDS-EC2), but consistent within the errors, as earlier found in other
studies (e.g. Shajib et al. 2020 and reference therein); 4) for KIDS-EC1 we
measure a quite strong external shear ($\gamma_{\rm ext}\sim 0.25$),
compatible with a group/cluster potential which is confirmed by the presence
of 12 more galaxies at a similar redshift as the lens within
$30^{\prime\prime}$ distance in the MUSE FOV (Napolitano et al. in prep.); 5)
both sources are aligned with the lens center of mass within 0.1′′, while the
stellar and mass centers are consistent within the errors. In Appedix we
discuss the impact of more general model density assumptions and show that
these do not impact the main conclusions of this study. From the lens model
parameters, we infer a projected mass within $R_{\rm E}$ of $\log M(R_{\rm
E})/M_{\sun}=11.28\pm 0.02$ for KIDS-EC1 and $\log M(R_{\rm
E})/M_{\sun}=11.42\pm 0.01$ for KIDS-EC2 (see also Table 2). Due to the
underlying assumption of a SIE mass distribution, we can easily derive the
mass inside $R_{\rm e}/2$ (i.e., $M(R_{\rm e}/2)$ in Table 2) to compare with
the dynamical mass by the velocity dispersion measurements inside the same
radius, derived above. We remark here that the adoption of $R_{\rm e}/2$ as
reference radius is consistent with previous strong lensing studies at the
same scale (e.g. Auger et al. 2010).

For the dynamical masses, we use the projected solution of the Jeans Equation
inside a circular aperture (e.g. Tortora et al. 2009), in order to take
correctly into account the light profile of the lenses (i.e. $n-$index)333This
is the most accurate way to determine the mass inside an aperture and avoid
assumptions about the virial estimates for a non-de Vaucouleurs profile
($n\neq 4$), see e.g. Cappellari et al. 2006. These projected masses, $M_{\rm
J}$ in Table 2, are fully consistent with the equivalent lensing-derived
masses, hence confirming the self-consistency of our mass estimates of the two
systems.

We finally estimate the total stellar mass, $M_{*},$ of the two lens systems
via SED fitting of the 9-band photometry, given in Table 1. We use the public
SED fitting–Code Investigating GALaxy Emission (Cigale v2018.0, Boquien et al.
2019. The $M_{*},$ are used, in combination with the lensing masses, to derive
the DM fraction inside $Re/2$, $f_{\rm DM}=1-M_{*}(R_{\rm e}/2)/M(R_{\rm
e}/2)$. We assume solar metallicity, while all other parameters, such as the
e-folding time, age of the main stellar population and internal extinction,
E(B-V), are free to vary. For the star formation rate (SFR), we adopt a
delayed star formation history, which allows us to efficiently model both
typical early-type and late-type galaxies (see Boquien et al. 2019). To double
check the results we also use another independent code, Le Phare (Ilbert et
al. 2006), with a similar set-up and found consistent results for all the
constrained parameters, within the errors. The stellar mass inside $R_{\rm
e}/2$ is derived assuming a constant stellar mass-to-light ratio for each lens
Sérsic light profile, whose parameters have been inferred from the lensing
model in Table 2. The final estimates of $M_{*}(R_{\rm e}/2)$ are reported in
the same Table, where we also list the final $f_{\rm DM}=0.60\pm 0.07$ for
KIDS-EC1 and $f_{\rm DM}=0.56\pm 0.04$ for KIDS-EC2, and report the difference
between the lensing and dynamical DM fractions ($\Delta f_{\rm DM}$).

The inferred $f_{\rm DM}$ are typical of DM dominated systems and consistent
with previous estimates based on lensing (Auger et al. 2010, Tortora et al.
2010, Schuldt et al. 2019) or dynamics of local (Tortora et al. 2009,
Cappellari et al. 2012, Tortora et al. 2012) or higher redshift galaxies
(Beifiori et al. 2014, Tortora et al. 2018) of similar stellar masses.

##  3 Characterization of the sources

A striking outcome of the lensing model is the degree of compactness of the
source galaxies in both ECs, that turned out to have effective radii smaller
than 0.1′′. The lensing model also provides an estimate of the $n-$indexes
($<1$ in both cases) and axis ratios ($\sim 0.9$ and $\sim 0.6$ respectively),
suggesting that the two sources might be disk dominated systems444The best-fit
$n-$indexes are fairly small if compared to typical disks: however we have
checked that fixing $n=1$, the other parameters change within the errors and
the reduced $\chi^{2}$ is worsened, hence demonstrating that $n-$indexes are
realistically $\;\hbox to0.0pt{\lower 2.5pt\hbox{$\sim$}\hss}\raise
1.5pt\hbox{$<$}\;1$ (see e.g. former findings on lensed quenching galaxies by
Geier et al. 2013).. Only space or adaptive-optics imaging will provide
accurate constraints for these parameters and confirm these findings, but
meanwhile, to better assess the reliability of the ground-based inferences, we
test the procedure using 20 mock ECs. We follow the same approach used to
simulate the lensing systems centered on randomly selected red luminous
galaxies used to train our CNN on KiDS ground-based images (see Li+20). In
particular we produce quad configurations using a source effective radius
varying between $R_{\rm e}=$[0.05′′, 0.15′′], $n-$index=[0.05,2], lens
effective radius $R_{\rm e}=$[1.0′′, 1.5′′] and Einstein radius, $R_{\rm
E}$=[1.0′′, 1.5′′]. The ranges adopted are meant to cover the parameter space
embracing the two crosses, and to demonstrate that the lens model tool can
recover the parameters correctly. In particular we have tested the case of
very small $R_{\rm e}$ and $n-$index of the source by modeling 4/20 mock
lenses with $R_{\rm e}<0.07^{\prime\prime}$ and $n-$index$<0.6$. The simulated
lenses are then convolved with a typical $r-$band PSF of KiDS observations and
noise is finally added to produce realistic KiDS-like r-band EC images. We
then run the lfit$\\_$gui using the same configuration file on these mock ECs
and derived the lensing parameters like for the real ECs. The derived source
$R_{\rm e}$ and $n-$index fall on the one-to-one relation with the input ones.
To quantify this we derive the following quantities: $\Delta R_{\rm e}/R_{\rm
e}=(R_{\rm e,in}-R_{\rm e,out})/R_{\rm e,in}=-0.01\pm 0.06$ and $\Delta
n/n=(n_{\rm in}-n_{\rm out})/n_{\rm in}=-0.01\pm 0.15$. The scatter is even
smaller ($\sim 0.03$ both in $R_{\rm e}$ and $n$) for the 4 most extreme
cases. Since they are both consistent with zero and the ($1\sigma$) scatter is
consistent with typical errors from the best-fit parameters in Table 2
($\;\hbox to0.0pt{\lower 2.5pt\hbox{$\sim$}\hss}\raise
1.5pt\hbox{$<$}\;20\%$), we are confident that the “compactness” of the
sources as well as its “disk-like” nature are real, although for the latter
there might be more freedom about the exact value of the $n-$index. We also
remark that sub-pixel sizes of the sources are refereed to the source plane,
where lfit$\\_$gui maps the source model with a spatial resolution 10$\times$
higher than the pixel scale of the lens plane, where the source images are
observed ($0.2^{\prime\prime}$ for KiDS images). Hence sizes smaller than a
single pixel are well within the reach of the tool we use. However, we have
found that 0.04′′ is the lowest limit for our groud-based observations. In
fact, mock lenses having sources with $R_{\rm e}$ in the range [0.01,0.04] are
recovered with too large uncertainties.

As described in §2.1, the source spectra are dominated by absorption lines,
which are typical of a relatively evolved population, and show only weak
emission lines (see Fig. 1), somehow at odds with the blue colors produced by
their continuum. In particular, KIDS-EC1 shows a low-SNR [OII] emission, while
KIDS-EC2 shows [OII] and H$\epsilon$, H$\delta$ and H$\gamma$ Balmer lines in
emission, although the latter are superimposed to absorption Balmer lines,
which makes any modeling of their profile very degenerate (e.g. continuum
level, velocity dispersion of the absorption and emission lines, age,
metallicity). Using the [OII] line, after some careful continuum subtraction,
we estimate a tentative SFR of $\sim 0.1M_{\sun}/$yr (using Eq. 3 in Kennicutt
1998). We stress that a SNR$\sim 6$ is not high enough for a robust estimate.

To better characterise the nature of the two sources, we thus exploit their
9-band photometry (Table 1) and run Cigale, as already done for the
deflectors.  In order to increase the SNR of the source SEDs, we average the
fluxes of the higher magnified lensed images (AB for KIDS-EC1 and AC for KIDS-
EC2) and obtain the de-lensed (i.e. using the $r-$band magnification,
$\mu_{\rm r}$, as in Table 2) SED in Fig. 3. We try two extreme metallicity
scenarios: a standard solar metallicity and a largely sub-solar one
($Z=0.0004$), being the latter suggested by the strength of the Balmer lines
compared with the Calcium H and K of KIDS-EC2 (for EC1 these lines are too
redshifted to fall in the MUSE wavelength range). For KIDS-EC1 we obtain a
better fit with solar metallicity (reduced $\chi^{2}=1.7$, corresponding to a
$\sim$10% significance for 8 degree-of-freedom), while for EC2 the fit is
better for the sub-solar metallicity ($\chi^{2}=0.7$, i.e. $\sim$ 70%
significance). The best fits, shown in Fig. 3, give for KIDS-EC1 an age of
$4.0\pm 0.1$ Gyr and a SFR of $9.1\pm 1.2$ $M_{\sun}/$yr, while for for KIDS-
EC2, age=$3.8\pm 0.4$ Gyr and SFR$=0.4\pm 0.1$ $M_{\sun}/$yr555For
completeness, the results obtained for KIDS-EC1 assuming a sub-solar
metallicity are age=$4.0\pm 1.0$ Gyr and SFR$=0.4\pm 0.1$ $M_{\sun}/$yr and
the results for KIDS-EC2 with solar metallicity are age=$1.2\pm 0.2$ Gyr and
SFR$=5.5\pm 0.3$ $M_{\sun}/$yr. Note that full spectro-photometric stellar
population analysis of the sources is beyond the purpose of this letter and
will be addressed in a separate paper.. Hence both systems have an old age for
their redshift, but are still forming stars at a low rate, although we see
only few emission lines in their spectra. In both cases the sources are quite
massive: KIDS-EC1 has a stellar mass of $\log
M_{*}/M_{\sun}=11.08^{+0.04}_{-0.05}$ and KIDS-EC2 of $\log
M_{*}/M_{\sun}=10.21^{+0.02}_{-0.02}$ , hence their specific star formation
rate (sSFR) are $\log{\rm sSFR}/$Gyr$=-1.12\pm 0.07$ and $\log{\rm
sSFR}/$Gyr$=-1.61\pm 0.03$ respectively, i.e. generally lower than the the
typical values expected for the main sequence (MS) of star-forming (SF)
galaxies at $z>1$ in the same mass range ($\log{\rm sSFR}_{\rm MS}/$Gyr$\sim
0.0\pm 0.4$, see e.g. Johnston et al. 2015).

Also for the lensed images we have performed a double check with Le Phare and
we have found slightly lower masses (0.1 and 0.2 dex for KIDS-EC1 and EC2,
respectively) but consistent or even lower $\log{\rm sSFR}/$Gyr ($-1.3$ and
$-1$ respectively), hence confirming the quenching status of the two systems.

Are these systems special? They are for their sizes, as they are both outliers
of the typical size–mass relation of SF systems by $\sim 0.5$dex in effective
radius. Allen et al. (2017), for instance, found a mean $\log R_{\rm
e}/$kpc$\sim 0.6$ and $0.5$ in the redshift bin $1<z<1.5$ for $\log
M_{*}/M_{\sun}\sim 11.1$ and $\sim 10.2$ respectively, while we have $\log
R_{\rm e}/$kpc$=-0.06$ for KIDS-EC1 and $-0.34$ for KIDS-EC2. This means that
they deviate significantly from normal galaxies at their redshift, while they
are closer to typical sizes of SF galaxies at $z~{}\hbox
to0.0pt{$>$\hss}{\lower 4.30554pt\hbox{$\sim$}}5$. On the other hand, they
also show very low sSFR, deviating from the MS by $\Delta_{\rm
MS}=\log$sSFR${}_{\rm MS}-\log$sSFR$\sim-1.1$ and $-0.4$. Simulations from
Tacchella et al. (2016) have shown that these low sSFR, together with sizes of
$\log R_{\rm e}/$kpc$<0$, are typical of a “post–blue nugget” (BNs, hereafter)
phase, i.e. systems having gone through compaction and entering their
quenching phase. Huertas-Company et al. (2018) have also found that most of
the massive, compact systems at z$>1$ ($\log M_{*}/M_{\sun}~{}\hbox
to0.0pt{$>$\hss}{\lower 4.30554pt\hbox{$\sim$}}10.3$) tend to be in such a
post–BN phase.

##  4 Conclusions and Perspectives

We have presented the confirmation and modelling of two Einstein Crosses found
in the KiDS footprint. The confirmation is based on MUSE spectroscopy,
unequivocally showing the strong-lensing nature of the systems. We have
detected the same spectral features in the four images of the two crosses and
inferred a redshift of $z_{\rm s}\sim 1.59$ for KIDS-EC1 and $z_{\rm s}\sim
1.10$ for KIDS-EC2, both higher that their respective deflector galaxies, i.e.
two old early-type galaxies at $z_{\rm l}\sim 0.38$ (KIDS J232940-340922) and
$z_{\rm l}\sim 0.24$ (KIDS J122456+005048).

The discovery is exceptional as we only inspected $\sim 1000$deg2 so far and
Einstein Crosses (ECs) are very rare phenomena. However, general predictions
on the number of expected quads (of which ECs are a special case) are based on
the assumption that these are generated by quasars (see e.g. Oguri & Marshall
2010). For the two ECs presented here, we have shown that they are produced by
ultra–compact, blue, quenching galaxies. In particular:

  1. 1.

to reproduce the cross configuration, the best lensing model (see Table 2),
assuming a singular isothermal ellipsoid (SIE) with external shear, predicts
that the sources have very compact sizes: $R_{\rm e}\sim 0.9$ kpc for KiDS-EC1
and $R_{\rm e}\sim 0.5$ kpc for KiDS-EC2, i.e. $>10\sigma$ off the typical
size-mass relation of normal SF galaxies at the same redshifts (see Allen et
al. 2017);

  2. 2.

the spectra of the sources show a dominance of absorption lines, typical of a
quite evolved stellar population, which has been confirmed by SED fitting of
the 9-band photometry performed on the average of the two highest magnified
images of each system. We inferred old ages (4 Gyr and 3.8 Gyr for KIDS-EC1
and KIDS-EC2, respectively) and moderate star formation ($9.1$ and
$0.4~{}M_{\sun}/$yr, respectively). However, the inferred stellar masses
($\log M_{*}/M_{\sun}\sim$11.08 and 10.21, respectively) imply very low
specific star-formation rates ($\log$sSFR$=-1.1$ and -1.6 Gyr-1,
respectively), typical of quenching galaxies at $z>1$ (see e.g. Tacchella et
al. 2016, Newman et al. 2018).

In §3, we have argued that the combination of an extremely compact size and
the low sSFR is compatible with the sources being two massive post–BNs. These
are compact massive galaxies having almost exhausted their star-forming phase
and currently undergoing quenching. As such, these systems are important to
understand the transformation of primordial disks into the compact cores (“red
nuggets”, see Dekel & Burkert 2014) of today’s large elliptical galaxies in
the first phase of their evolution, before they enter their subsequent merging
phase (e.g. Oser et al. 2012). These systems have been predicted to be very
numerous in simulations (e.g. Zolotov et al. 2015, Tacchella et al. 2016) and
their census at $z>1$ has just started, including the confirmation of their
abundance and physical properties (see e.g. Huertas-Company et al. 2018).
Having observed two of such systems in a peculiar lensing configuration, might
suggest that, indeed, they are not uncommon at $z>1$. To make a rough estimate
of the expected numbers of such EC events from compact post–BN systems, we can
use predictions based on a size-luminosity relation compatible with high$-z$
studies calibrated over compact star-forming galaxies from Collett (2015),
reproducing typical ground based surveys with seeing $\leq
1^{\prime\prime}$666We have used the pre-compiled predictions for the Dark
Energy Survey (DES), with seeing 1′′ in github.com/tcollett/LensPop. These
simulations are consistent with KiDS in terms of depth ($r\sim$25, 5$\sigma$
within 2′′), hence giving access to a similar lens luminosity distribution,
and only slightly worse in image quality ($\sim 1^{\prime\prime}$ for DES vs
0.7′′ for KiDS and a pixel size 0.26′′ vs 0.21′′). Furthermore the adoption of
a slightly larger seeing allows us to better account for the selection
function introduced by the visual inspection, as the human eye tends to give a
lower grade to arcs or multiple images that are too close to the lens center
as they look diluted/confused in the lens starlight.. We estimate $\sim 1$ EC
over 1000deg2 (assuming $<0.1^{\prime\prime}$ alignment between source and
deflector), generated by a slightly flattened mass distribution like the one
estimated for our ECs ($q=$[0.7, 0.9]) by a compact source (size $\;\hbox
to0.0pt{\lower 2.5pt\hbox{$\sim$}\hss}\raise 1.5pt\hbox{$<$}\;1$ kpc) at
$z<1.7$, with sufficient SNR ($>10$) to be identified as a sure lens around
bright deflectors ($r<21$, according to our selection in Li+20). This rough
prediction can be likely an upper limit because we are assuming that all
sources aligned within 0.1′′ produce an EC, but using mock lenses we have
checked that most of them should look as such (see §3). The two ECs reported
here are possibly slightly overabundant with respect to the expectation for
standard star-forming systems. This is a first interesting indication that
this population of post–BNs might be particularly abundant at $z>1$. According
to the same predictions, releasing any limitation on the brightness of the
lens (i.e. $r<25$), the expected EC/quad configurations are $\sim$ half a
dozen for sources at $z<4$ (see also §1). Hence we expect to confirm other
systems in KiDS, in the future. More importantly we can expect that for next
generation surveys like LSST, EUCLID or CSST, we can discover $4000$ to $7000$
similar quad-like configurations. This will represent a unique opportunity to
perform systematic studies of this population of compact systems in great
details using lensing as a “gravitational telescope”, which will be
prohibitive to observe without lensing magnification (see e.g. Toft et al.
2017) before extremely large telescopes will be online.

## Acknowledgements

We thank the anonymous referee for the constructive reports that helped us to
strengthen the results of our paper. NRN acknowledges financial support from
the “One hundred top talent program of Sun Yat-sen University” grant N.
71000-18841229. RL acknowledge support from China Postdoctoral Science
Foundation 2020M672935 and Guangdong Basic and Applied Basic Research
Foundation 2019A1515110286. RL also acknowledge the Fundamental Research Funds
for the Central Universities, Sun Yat-sen University 71000-31610034. CS is
supported by a Hintze Fellowship at the Oxford Centre for Astrophysical
Surveys. CT acknowledges funding from the INAF PRIN-SKA 2017 program
1.05.01.88.04. GD acknowledges support from CONICYT project Basal AFB-170002.
KK acknowledges support by the Alexander von Humboldt Foundation. CH
acknowledges support from the European Research Council under grant number
647112, and support from the Max Planck Society and the Alexander von Humboldt
Foundation. HH is supported by a Heisenberg grant of the Deutsche
Forschungsgemeinschaft (Hi 1495/5-1) and ERC Consolidator Grant (No. 770935).
MB is supported by the Polish Ministry of Science and Higher Education through
grant DIR/WK/2018/12, and by the Polish National Science Center through grants
no. 2018/30/E/ST9/00698 and 2018/31/G/ST9/03388. We thank Yiping Shu for
providing the lens fitting code lfit$\\_$gui, and for technical support. We
also thank Felipe Barrientos and Nicolas Tejos for useful discussions.

## References

  * Allen et al. (2017) Allen, R. J., Kacprzak, G. G., Glazebrook, K., et al. 2017, ApJ, 834, L11
  * Auger et al. (2010) Auger, M. W., Treu, T., Bolton, A. S., et al. 2010, ApJ, 724, 511
  * Beifiori et al. (2014) Beifiori, A., Thomas, D., Maraston, C., et al. 2014, ApJ, 789, 92
  * Bettoni et al. (2019) Bettoni, D., Falomo, R., Scarpa, R., et al. 2019, ApJ, 873, L14
  * Bolton et al. (2006) Bolton, A. S., Moustakas, L. A., Stern, D., et al. 2006, ApJ, 646, L45
  * Boquien et al. (2019) Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000
  * Burchett et al. (2020) Burchett, J. N., Rubin, K. H. R., Prochaska, J. X., et al. 2020, arXiv:2005.03017
  * Cappellari et al. (2006) Cappellari, M., Bacon, R., Bureau, M., et al. 2006, MNRAS, 366, 1126
  * Cappellari et al. (2012) Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2012, Nature, 484, 485
  * Cappellari (2017) Cappellari, M. 2017, MNRAS, 466, 798
  * Chabrier (2003) Chabrier, G. 2003, ApJ, 586, L133
  * Collett (2015) Collett, T. E. 2015, ApJ, 811, 20
  * de Jong et al. (2015) de Jong, J. T. A., Verdoes Kleijn, G. A., Boxhoorn, D. R., et al. 2015, A&A, 582, A62
  * Dekel & Burkert (2014) Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870
  * Edge et al. (2013) Edge, A., Sutherland, W., Kuijken, K., et al. 2013, The Messenger, 154, 32
  * Geier et al. (2013) Geier, S., Richard, J., Man, A. W. S., et al. 2013, ApJ, 777, 87
  * Huertas-Company et al. (2018) Huertas-Company, M., Primack, J. R., Dekel, A., et al. 2018, ApJ, 858, 114
  * Johnston et al. (2015) Johnston, R., Vaccari, M., Jarvis, M., et al. 2015, MNRAS, 453, 2540
  * Kelly et al. (2015) Kelly, P. L., Rodney, S. A., Treu, T., et al. 2015, Science, 347, 1123
  * Kennicutt (1998) Kennicutt, R. C. 1998, ARA&A, 36, 189
  * Koopmans et al. (2006) Koopmans, L. V. E., Treu, T., Bolton, A. S., et al. 2006, ApJ, 649, 599
  * Kormann et al. (1994) Kormann, R., Schneider, P., & Bartelmann, M. 1994, A&A, 284, 285
  * Kuijken et al. (2019) Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2
  * Ilbert et al. (2006) Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841
  * Li et al. (2020) Li, R., Napolitano, N. R., Tortora, C., et al. 2020, ApJ, 899, 30 (Li+20)
  * Magain et al. (1988) Magain, P., Surdej, J., Swings, J.-P., et al. 1988, Nature, 334, 325
  * Moré (1978) Moré, J. J. 1978, Lecture Notes in Mathematics, Berlin Springer Verlag, 105
  * Muzzin et al. (2012) Muzzin, A., Labbé, I., Franx, M., et al. 2012, ApJ, 761, 142
  * Newman et al. (2018) Newman, A. B., Belli, S., Ellis, R. S., et al. 2018, ApJ, 862, 125
  * Nightingale et al. (2019) Nightingale, J. W., Massey, R. J., Harvey, D. R., et al. 2019, MNRAS, 489, 2049
  * Oser et al. (2012) Oser, L., Naab, T., Ostriker, J. P., et al. 2012, ApJ, 744, 63
  * Petrillo et al. (2017) Petrillo, C. E., Tortora, C., Chatterjee, S., et al. 2017, MNRAS, 472, 1129
  * Petrillo et al. (2019a) Petrillo, C. E., Tortora, C., Vernardos, G., et al. 2019, MNRAS, 484, 3879
  * Petrillo et al. (2019b) Petrillo, C. E., Tortora, C., Chatterjee, S., et al. 2019, MNRAS, 482, 807
  * Oguri & Marshall (2010) Oguri, M., & Marshall, P. J. 2010, MNRAS, 405, 2579
  * Ostrovski et al. (2018) Ostrovski, F., Lemon, C. A., Auger, M. W., et al. 2018, MNRAS, 473, L116
  * Roy et al. (2018) Roy, N., Napolitano, N. R., La Barbera, F., et al. 2018, MNRAS, 480, 1057
  * Rubin et al. (2014) Rubin, K. H. R., Prochaska, J. X., Koo, D. C., et al. 2014, ApJ, 794, 156
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Shu et al. (2016) Shu, Y., Bolton, A. S., Mao, S., et al. 2016, ApJ, 833, 264
  * Schuldt et al. (2019) Schuldt, S., Chirivì, G., Suyu, S. H., et al. 2019, A&A, 631, A40
  * Shajib et al. (2020) Shajib, A. J., Treu, T., Birrer, S., et al. 2020, arXiv:2008.11724
  * Spiniello et al. (2011) Spiniello, C., Koopmans, L. V. E., Trager, S. C., et al. 2011, MNRAS, 417, 3000
  * Spiniello et al. (2018) Spiniello, C., Agnello, A., Napolitano, N. R., et al. 2018, MNRAS, 480, 1163
  * Spiniello et al. (2019) Spiniello, C., Agnello, A., Sergeyev, A. V., et al. 2019, MNRAS, 483, 3888
  * Suyu et al. (2013) Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70
  * Tacchella et al. (2016) Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 457, 2790
  * Tessore et al. (2016) Tessore, N., Bellagamba, F., & Metcalf, R. B. 2016, MNRAS, 463, 3115
  * Toft et al. (2017) Toft, S., Zabl, J., Richard, J., et al. 2017, Nature, 546, 510
  * Tortora et al. (2009) Tortora, C., Napolitano, N. R., Romanowsky, A. J., et al. 2009, MNRAS, 396, 1132
  * Tortora et al. (2010) Tortora, C., Napolitano, N. R., Romanowsky, A. J., et al. 2010, ApJ, 721, L1
  * Tortora et al. (2012) Tortora, C., La Barbera, F., Napolitano, N. R., et al. 2012, MNRAS, 425, 577
  * Tortora et al. (2018) Tortora, C., Napolitano, N. R., Roy, N., et al. 2018, MNRAS, 473, 969
  * Treu et al. (2010) Treu, T., Auger, M. W., Koopmans, L. V. E., et al. 2010, ApJ, 709, 1195
  * Zolotov et al. (2015) Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327

## Impact of the SIE assumption and the fitting tool

The mismatch between the stellar velocity dispersion measured from the
spectrum and the lensing dispersion of KIDS-EC1 found in §2.2 suggests that
the assumption of a SIE mass density might be inappropriate. To investigate
the impact of this assumption and check the reliability of the main parameters
derived by lfit_gui and reported in Table 2, we use a different lens model
tool, LENSED (Tessore et al. 2016). We stress that a full comparison of
different fitting techniques is beyond the purpose of this paper, however we
can use the relevant outputs of an independent tool to validate the main
results of this paper. First, as LENSED allows a free slope for the mass
density profile (i.e. an Elliptic Power Law, EPL, model), we can check how
this deviates from SIE (3D slope$=2$). For KIDS-EC1 we find the best-fit 3D
slope$=-1.75\pm 0.06$, which is consistent with the Auger et al. (2010)
formula ($\gamma^{\prime}-2=2.67(f_{\rm SIE}-1)+0.20$, where $\gamma^{\prime}$
is the 3D slope and $f_{\rm SIE}=\sigma_{\rm e/2}/\sigma_{\rm SIE}$) for the
mismatch between stellar velocity dispersion measured from the spectrum and
the lensing dispersion for the SIE model. For EC2 we have a 3D slope$=2.04\pm
0.01$, which is instead very close to the 3D slope predicted for a SIE. Due to
the different slopes, we derive different 2D mass extrapolation at $R_{\rm
e}/2$: for KIDS-EC1 we have $\log M(R_{\rm e}/2)/M_{\sun}=10.99\pm 0.07$, for
KIDS-EC2 we have $\log M(R_{\rm e}/2)/M_{\sun}=11.20\pm 0.02$. The inferred
dark matter fractions are $f_{\rm DM}=0.54\pm 0.11$ for KIDS-EC1 and $f_{\rm
DM}=0.58\pm 0.05$ for KIDS-EC2, i.e. fully consistent with the SIE results
within $1\sigma$. This shows that the inner slope does not impact
significantly our inference on the DM content of the galaxies, as well as the
overall conclusions related to the lensing parameters in Table 2, hence the
use of the SIE in the rest of the paper is fairly justified.

Another central result of the paper is the compactness of the source, which,
as discussed in §3, is challenging to assess from ground based imaging and
might depend, e.g., on the way a given tool samples the model at the sub-pixel
scale, performs the convolution with the PSF and re-bins the model. We have
double checked the results of lfit_gui against LENSED and confirmed both
compactness and the pseudo-exponential light profiles of the sources, even
though the central values are consistent only within $\sim 2\sigma$. In fact,
from LENSED we found that the source sizes of KIDS-EC1 is $R_{\rm
e}=0.05^{\prime\prime}\pm 0.01^{\prime\prime}$ while its $n-$index is
$n=1.2\pm 0.3$; for KIDS-EC2, we found instead $R_{\rm
e}=0.03^{\prime\prime}\pm 0.01^{\prime\prime}$ and $n=1.58\pm 0.07$. Finally,
source axis ratios and magnitudes are fairly similar to lfit_gui inferences:
for KiDS-EC1, LENSED gives ($b/a,~{}r$) $=$(0.24, 23.54), for KiDS-EC2 (0.64,
23.80). As comparison, lfit_gui gives (0.5, 23.78), (0.7, 23.69) for KiDS-EC1
and KiDS-EC2, respectively.

Table 1: Einstein Cross Optical-NIR photometry

ID |  $\Delta$RA |  $\Delta$DEC | $u$ | $g$ | $r$ | $i$ | $Z$ | $Y$ | $J$ | $H$ | $K_{\rm s}$
---|---|---|---|---|---|---|---|---|---|---|---
| (arcsec) | (arcsec) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag)
KIDS-EC1: KIDS J232940-340922 RA$=352.417753$ DEC$=-34.156375$
G | $0.00\pm 0.03$ | $0.00\pm 0.03$ | 21.44 | 20.69 | 19.77 | 19.29 | 19.05 | 18.48 | 17.90 | 17.19 | 16.48
|  |  | (0.22) | (0.09) | (0.04) | (0.04) | (0.12) | (0.06) | (0.08) | (0.08) | (0.04)
A | $-1.04\pm 0.06$ | $-0.75\pm 0.05$ | 22.57 | 21.97 | 22.29 | 22.27 | 21.16 | 21.14 | 20.46 | 19.70 | 19.39
|  |  | (0.17) | (0.15) | (0.11) | (0.17) | (0.11) | (0.14) | (0.12) | (0.12) | (0.12)
B | $0.87\pm 0.08$ | $0.99\pm 0.05$ | 22.65 | 21.98 | 22.36 | 22.10 | 21.32 | 21.29 | 20.24 | 19.60 | 19.44
|  |  | (0.19) | (0.11) | (0.12) | (0.15) | (0.13) | (0.17) | (0.10) | (0.10) | (0.12)
C | $-0.58\pm 0.09$ | $0.67\pm 0.07$ | 24.06 | 22.24 | 22.62 | 22.23 | 21.31 | 21.18 | 20.38 | 20.14 | 20.32
|  |  | (0.79) | (0.17) | (0.29) | (0.37) | (0.24) | (0.21) | (0.16) | (0.23) | (0.39)
D | $0.61\pm 0.19$ | $-0.72\pm 0.16$ | 25.44 | 22.80 | 23.14 | 23.24 | 21.86 | 21.81 | 21.32 | 20.23 | 19.57
|  |  | (1.65) | (0.39) | (0.38) | (0.51) | (0.39) | (0.71) | (0.42) | (0.43) | (0.18)
KIDS-EC2: KIDS J122456+005048 RA$=186.233401$ DEC$=+0.846682$
G | $0.00\pm 0.01$ | $0.00\pm 0.01$ | 20.94 | 19.30 | 18.10 | 17.48 | 17.20 | 17.04 | 16.67 | 16.24 | 15.98
|  |  | (0.09) | (0.03) | (0.03) | (0.04) | (0.03) | (0.05) | (0.06) | (0.02) | (0.38)
A | $-1.24\pm 0.05$ | $0.78\pm 0.05$ | 23.13 | 22.52 | 21.90 | 21.57 | 20.74 | 20.51 | 20.34 | 20.33 | 19.82
|  |  | (0.14) | (0.16) | (0.17) | (0.13) | (0.08) | (0.09) | (0.11) | (0.08) | (0.42)
B | $0.50\pm 0.07$ | $1.29\pm 0.08$ | 23.25 | 23.14 | 22.01 | 21.79 | 21.05 | 20.97 | 20.66 | 20.56 | 19.87
|  |  | (0.11) | (0.33) | (0.21) | (0.21) | (0.09) | (0.11) | (0.12) | (0.09) | (0.42)
C | $1.27\pm 0.06$ | $-0.96\pm 0.06$ | 23.44 | 22.63 | 21.97 | 21.65 | 20.73 | 20.62 | 20.43 | 20.08 | 19.96
|  |  | (0.15) | (0.16) | (0.15) | (0.13) | (0.08) | (0.11) | (0.11) | (0.08) | (0.44)
D | $-0.73\pm 0.09$ | $-1.14\pm 0.09$ | 23.36 | 22.77 | 22.11 | 21.93 | 21.22 | 21.12 | 20.94 | 20.69 | 20.34
|  |  | (0.18) | (0.29) | (0.22) | (0.27) | (0.15) | (0.16) | (0.19) | (0.10) | (0.46)

Note. — Objects coordinates are in degree, errors on magnitudes in different
bands are given in brackets. Magnitudes of the lensed images (A, B, C, D) are
total magnitudes obtained by Gaussian fit after the central galaxy G has been
removed. Galaxy magnitude are total Sérsic magnitudes.

Table 2: Summary of the Einstein cross main parameters

Parameter | KIDS-EC1 | KIDS-EC2
---|---|---
MUSE spectroscopy
$z_{\rm l}$ | $0.3810\pm 0.001$ | $0.2372\pm 0.0005$
$z_{\rm s}$ | $1.59\pm 0.01$ | $1.10\pm 0.01$
$\sigma_{\rm R_{\rm e}/2}$ (km s-1) | $192\pm 4$ | $248\pm 5$
Lensing model
$\theta_{\rm E}$ (arcsec) | $0.99\pm 0.02$ | $1.42\pm 0.01$
$R_{\rm E}$ (kpc) | $5.2\pm 0.1$ | $5.4\pm 0.1$
$\sigma_{\rm E}$ (km s-1) | $226\pm 8$ | $260\pm 5$
lens $b/a$ | $0.91\pm 0.07$ | $0.68\pm 0.07$
lens $PA$ | $1\pm 56$ | $148\pm 2$
$\mu_{r}$ (A, B, C, D) | (3.6, 3.4, 2.7, 1.7) | (5.6, 5.1, 5.3, 4.6)
lens $PA$ | $1\pm 56$ | $148\pm 2$
star $n-$index | $2.8\pm 0.3$ | $3.43\pm 0.08$
star $b/a$ | $0.89\pm 0.03$ | $0.59\pm 0.01$
star $PA$ | $107\pm 9$ | $146\pm 1$
star $R_{\rm e}$ (arcsec) | $1.12\pm 0.10$ | $1.63\pm 0.05$
star $R_{\rm e}$ (kpc) | $5.8\pm 0.5$ | $6.2\pm 0.2$
$R_{\rm E}/R_{\rm e}$ | $0.90\pm 0.08$ | $0.87\pm 0.03$
star $\Delta$RA | $-0.02\pm 0.02$ | $0.01\pm 0.01$
star $\Delta$DEC | $0.02\pm 0.02$ | $0.01\pm 0.01$
source $n-$index | $0.4\pm 0.1$ | $0.2\pm 0.1$
source $b/a$ | $0.5\pm 0.1$ | $0.7\pm 0.1$
source $R_{\rm e}$ (arcsec) | $0.10\pm 0.02$ | $0.06\pm 0.01$
source $R_{\rm e}$ (kpc) | $0.87\pm 0.14$ | $0.46\pm 0.04$
source $\Delta$RA | $-0.03\pm 0.02$ | $0.04\pm 0.01$
source $\Delta$DEC | $-0.07\pm 0.02$ | $0.05\pm 0.01$
external shear $\gamma_{\rm ex}$ | $0.25\pm 0.03$ | $0.01\pm 0.01$
Mass estimates
$\log M(R_{\rm E})/M_{\sun}$ | $11.28\pm 0.02$ | $11.42\pm 0.01$
$\log M(R_{\rm e}/2)/M_{\sun}$ | $11.04\pm 0.03$ | $11.18\pm 0.01$
$\log M_{\rm J}(R_{\rm e}/2)/M_{\sun}$ | $10.95\pm 0.05$ | $11.21\pm 0.02$
$\log M_{*}/M_{\sun}$ | $11.18\pm 0.04$ | $11.33\pm 0.03$
$\log M_{*}/M_{\sun}(R_{\rm e}/2)$ | $10.65\pm 0.04$ | $10.81\pm 0.03$
$f_{\rm DM}(R_{\rm e}/2)$ | $0.60\pm 0.07$ | $0.56\pm 0.04$
$\Delta f_{\rm DM}=f_{\rm DM,len}-f_{\rm DM,J}$ | $0.09\pm 0.10$ | $-0.02\pm 0.07$

Note. — Einstein cross parameters. Muse spectroscopy: redshift of the lens and
source and velocity dispersion of the lens calculated at $R_{\rm e}/2$.
Lensing models: Einstein radius in arcsec ($\theta_{\rm E}$) and kpc ($R_{\rm
E}$) and the model SIE velocity dispersion, $\sigma_{\rm E}$, followed by
self-explaining parameters related to the total mass (labelled by ‘lens’),
stellar mass parameters (labelled by ‘star’) and source parameters (labelled
by ‘source’) light distribution. Offsets ($\Delta$RA, $\Delta$DEC) are
calculated with respect to the lens center. Mass estimates: summary of the
mass estimates from lensing model ($M$), Jeans model ($M_{J}$), and stellar
population $M_{*}$ (see text for details) together with the projected DM
fractions ($f_{\rm DM}$) and difference between lens ($f_{\rm DM,len}$) and
Jeans ($f_{\rm DM,J}$) analyses.

Figure 1: Detection and confirmation of the Einstein crosses: KIDS
J232940-34092 (top), KIDS J122456+005048 (bottom). Left column: for each cross
we show the KiDS ($gri$) color image (top) and the MUSE white light image
(bottom). Right column: for each cross we plot the MUSE spectrum of the
deflector, in black, with the best fit model for the velocity dispersion
estimate with pPXF, in red, (top) and of the 4 individual images, in blue
tones, and the mean spectrum, in black (bottom). Over-plotted on all spectra,
are the main absorption lines (in red) and emission lines (in blue), shifted
to the estimated redshift of the represented object (see Table 2). Figure 2:
Lensing models for KIDS J232940-34092 (top), KIDS J122456+005048 (bottom).
From left to right we show the $r-$band KiDS image used for the model (A), the
foreground light-subtracted image (B), the Einstein cross images with the
foreground light subtracted (C), the reconstructed Einstein crosses, the
residual image (E=A-B-D), and the reconstructed background source (F). We use
a singular isothermal ellipsoid to model the deflector mass model and a Sérsic
model for foreground and source light. For KIDS J122456+005048, to model the
foreground light, we also account for the light of the two nearby galaxies
marked as G1 and G2 (see panel B). They are also visible in the KiDS and MUSE
images (see Fig. 1) and found to have a similar redshift as the lens. Figure
3: SED fitting of the mean optical+NIR photometry obtained from averaging the
de-lensed fluxes of the most magnified images in Table 1. KIDS J232940-34092
photometry (black points with error bars) is plotted together with the best
fit template model (gray solid line), corresponding to the parameters reported
in the figure (see text for more details) and solar metallicity. Dark blue
points show the corresponding photometry from the model used to fit the
observations. KIDS J122456+005048 photometry (red points with errorbars), is
plotted against the best fit model (light red line) and integrated photometry
(yellow points), as well as model parameters for the sub-solar metallicity
case (see text for details). Data from KIDS J232940-34092 have been shifted
+0.005 upward for clarity.