# Predicting Text Readability from Scrolling Interactions

Sian Gooding
Dept of Computer Science and Technology
University of Cambridge
<EMAIL_ADDRESS>
Yevgeni Berzak
MIT BCS
<EMAIL_ADDRESS>
Tony Mak
Google
<EMAIL_ADDRESS>
Matt Sharifi
Google
<EMAIL_ADDRESS>

###### Abstract

Judging the readability of text has many important applications, for instance
when performing text simplification or when sourcing reading material for
language learners. In this paper, we present a $518$ participant study which
investigates how scrolling behaviour relates to the readability of English
texts. We make our dataset publicly available and show that (1) there are
statistically significant differences in the way readers interact with text
depending on the text level, (2) such measures can be used to predict the
readability of text, and (3) the background of a reader impacts their reading
interactions and the factors contributing to text difficulty.111Our dataset is
available at https://github.com/siangooding/readability_scroll

##  1 Introduction

There are multiple attributes of a written text that impact how difficult it
is to read. This concept is formally known as the readability of text, where
readability is defined as the sum total of elements within a text that impact
a reader’s understanding, reading speed and level of interest Dale and Chall
(1949). Many factors can influence the readability of text, such as the
lexical and syntactic complexity, level of conceptual difficulty and style of
writing Xia et al. (2016). For instance, Figure 1 presents an example sentence
(b) that has been rewritten to a more readable format shown in (a).222Example
from the OneStopEnglish dataset by Vajjala and Lučić (2018).

| (a) A new scientific study says that global
---
warming might make temperatures rise
more than people think.
(b) Temperature rises resulting from unchecked climate change will be at the
severe end of those projected, according to a new scientific study.
Figure 1: Example of two sentences (a) and (b) with differing levels of
readability. Sentence (b) has a more sophisticated syntactic structure and
advanced vocabulary than (a).

Automatically measuring the readability of text has many useful applications.
For example, it is used in text simplification Aluisio et al. (2010), when
sourcing reading material for language learners Xia et al. (2016), in ranking
search engine content for dyslexic users Morris et al. (2018) and to ensure
critical consumer information is delivered at an appropriate level Zou et al.
(2019). However, current approaches for measuring readability rely exclusively
on linguistic features which do not account for the subjective needs of
readers. As a result, traditional readability formulas perform poorly in
modeling adult judgements of textual complexity Crossley et al. (2017).

Furthermore, traditional readability approaches do not work well for online
content Collins-Thompson (2013). This is due to systems being highly sensitive
to noise; requiring well formed sentences and performing poorly on short
passages Collins-Thompson (2014). Linguistically driven techniques are also
language specific and require sophisticated models to extract features. For
low-resource languages, such tools and models may be unavailable Agić et al.
(2016).

Achieving a more inclusive assessment of text readability requires collecting
subjective measures on multiple levels (e.g., reading interactions and
comprehension-based questions) across different genres and populations.
Therefore, alternative approaches must be examined to further the
understanding of text difficulty and the multiple factors that influence it.
In this paper, we introduce such an alternative approach, by using
crowdsourcing to collect scroll-based interactions whilst participants read
texts at differing levels. We present a $518$ participant study, recording the
reading interactions of individuals and obtain statistically significant
results showing that reading interactions differ dependent on the textual
complexity.

In the following sections, we discuss relevant previous work, outline our data
collection and report results from our study. As a preliminary evaluation, we
report the statistical significance of reading interactions based on text
difficulty, and integrate scroll features into a readability classifier.
Finally, we investigate group differences in readability based on the first
language of participants. To conclude, in this paper we make the following
contributions:

  * •

We release a dataset containing the reading interactions and comprehension
scores of $518$ participants based on the OneStopEnglishQA dataset Berzak et
al. (2020).

  * •

We show that scroll-based reading interactions can be used to predict the
readability of text as well as improve current state-of-the-art approaches.

  * •

We investigate how first language impacts reading interactions and emphasise
that the factors contributing to text difficulty vary depending on the target
audience.

##  2 Related Work

### Text Readability

Research on assessing the readability of English texts has spanned several
decades. The earliest works focused on the construction of readability
formulae and metrics Chall (1958); Klare et al. (1963); Zakaluk and Samuels
(1988). These measures rely on shallow textual characteristics such as the
number of sentences, average length of sentences and average word length. For
example, one of the most well known readability scores is the Flesch-Kincaid
score Kincaid et al. (1975). This score takes the average number of words per
sentence as well as the average number of syllables as a proxy for syntactic
and semantic difficulty. However, there are multiple limitations to
readability formulae Collins-Thompson (2014). Firstly, these formulae are
based on surface characteristics of text, and ignore deeper levels of text
processing known to be important factors of readability. Furthermore, these
metrics typically assume the text will contain no noise. As a result of this,
a number of studies have demonstrated that these metrics are unreliable for
web-based content Si and Callan (2001); Collins-Thompson and Callan (2004);
Feng et al. (2009).

Due to the limitations of traditional formulae, research in readability
assessment subsequently shifted towards machine learning (ML) based
techniques. These approaches combine a far richer variety of linguistic
features using ML classification algorithms and result in far better
performance François and Miltsakaki (2012). Early work on statistical
readability assessment demonstrated the improvement of these approaches for
readability prediction Si and Callan (2001); Collins-Thompson and Callan
(2004). Subsequent work focused on the addition of appropriate features, for
instance, lexical and grammatical Heilman et al. (2007, 2008) as well as
discourse-based Pitler and Nenkova (2008); Feng et al. (2010); Graesser et al.
(2011). Systems specifically designed for non-native audiences have also been
developed Feng et al. (2010); Vajjala and Meurers (2012); Xia (2019).
Additionally, the use of eye-tracking data from both language learners and
native speakers has been shown to improve readability assessment models
González-Garduno and Søgaard (2018). However, many of the aforementioned
models require extensive feature engineering, which not only rely on well-
formed content, but also depend on multiple resources such as parsers and
word-lists.

### Implicit Feedback Techniques

Reading on a device, such as a tablet or phone, has predominantly taken the
place of traditional formats Li et al. (2019). Such devices allow access to
implicit user feedback by measuring how a user interacts with the text they
read. A key advantage of implicit feedback techniques is that they can
unobtrusively obtain information by measuring user interactions with a system
Kelly and Teevan (2003). For instance, Claypool et al. (2001) designed a study
to capture mouse and keyboard interactions as implicit measures of interest.
They measure the correlation of these implicit features with user ratings and
find the time spent and amount of scrolling had a positive correlation with
interest. Chen et al. (2021) analyse the factors predictive of English
language typing times to investigate effects of linguistic context on language
production. Implicit feedback techniques have additionally been shown as
useful in information retrieval Golovchinsky et al. (1999), ranking summaries
White et al. (2002) and identifying user preferences Kelly and Teevan (2003).

In our study, we measure implicit feedback from participant interactions
whilst reading. The goal is to produce a set of readability features, using
aggregate scroll interactions, that are robust to noise and do not require
extensive feature engineering or external resources.

Figure 2: Example of the reading interface presented to participants.

Elementary
---
| The brown beer bottle was in the water for 101 years. A
---
fisherman found it in the Baltic Sea off the Northern city
of Kiel. Holger von Neuhoff a curator at the museum
said this bottled message was the oldest he had ever seen.
Advanced
The brown beer bottle, which had been in the water for 101 years, was found in
the catch of Konrad Fischer, a a fisherman, who had been out in the Baltic Sea
off of the northern city of Kiel. Holger von Neuhoff, curator for ocean and
science at the museum, said that this bottled message was the oldest he had
come across.
| Q: How was the bottle discovered?
---
a: A fisherman discovered it
b: It washed up on the shore of Kiel
c: Holger von Neuhoff discovered it
d: A boy discovered it while playing in the sea
Figure 3: Example paragraph from elementary and advanced texts. In the figure,
answer spans are highlighted in green and distractor spans in pink.

##  3 Data Collection

We design a controlled reading task and record user reading interactions. We
collect responses from participants in the United States and India. This
decision was made based on our pilot studies which showed participants from
India were much more likely to know English as a second language. The
experiment was designed in line with common practises in scroll interaction
research Hinckley et al. (2002). This study represents best efforts in the
collection of a novel dataset of reading interactions for text at different
levels of readability.

### Sampling and Participants

Participants in the study were recruited via a crowdsourcing platform and were
based in either India or the US. A range of background information pertinent
to language proficiency was collected using a demographic questionnaire which
is available on our dataset repository. This included the self-reported
English proficiency, native language, hours spent reading per week and highest
level of formal education. The design of the questionnaire was informed
following guidelines on judging reading ability Acheson et al. (2008). English
was the first language of $66$% of participants and a total of $46$ distinct
first languages were recorded. The second most frequent first language was
Tamil, accounting for $19$% of participants.

### Materials

Investigating how readability impacts user interactions required texts at
different levels of difficulty. For this, we use the OneStopEnglishQA dataset
Berzak et al. (2020) which has been designed for the evaluation and analysis
of reading comprehension in English. The dataset contains $30$ articles from a
prior dataset collected for readability assessment and automatic text
simplification: OneStopEnglish Vajjala and Lučić (2018). Each article has been
rewritten by teachers to suit three levels of adult ESL learners (elementary,
intermediate, and advanced) and are originally from the Guardian
newspaper.333https://www.theguardian.com/uk OneStopEnglishQA additionally
contains high quality text comprehension questions that are explicitly linked
to textual spans. In our study, we use only the advanced and elementary texts
and ask the same questions per article independent of the level. Figure 3
illustrates an example of a comprehension question and the relevant textual
spans for both an advanced and elementary paragraph. For each text, we present
three comprehension questions that have been strategically chosen to
necessitate scrolling. The texts presented to participants are chosen
randomly, however a participant will never be shown the same article at
differing levels.

### Task and Procedure

The task required participants to read articles and then answer three
comprehension questions. Initially, participants were shown instructions as
well as a placeholder text with questions before confirming that they were
ready to begin. Participants were presented with two texts, one at an advanced
level and one at an elementary level. Texts were shown one at a time in a
random order. Participants were immediately given the comprehension questions
after reading the text. The instructions stated that participants must read
the article and then answer the comprehension questions. However, the text was
still accessible during question answering and we record these interactions
separately to reading. Additionally, participants were informed that they
would be awarded a financial bonus on top of their base rate of pay for each
question answered correctly. The decision to award a bonus was to encourage
participants to read the text carefully. When presented with an article,
participants were not able to progress until at least $90$ seconds had elapsed
to encourage engagement and prevent immediate skipping. Once the participant
had read the two texts and answered all questions, a demographic questionnaire
was presented.

### Implementation

We used the Qualtrics platform to create the survey and added a custom front
end implementation using HTML, CSS and JavaScript. The experiment interface is
illustrated in Figure 3 and was displayed via a browser window. The article
text was presented to participants in a restricted window of size $1080$ by
$1920$ CSS px, a density-independent measure, based on the dimensions of an
average Android device. The font used in the study was sans-serif size $18$pt.

We log an event whenever the participant scrolls on the text. Events are
logged every $100ms$ and include a timestamp indicating the elapsed time as
well as the scroll y-axis offset in pixels (the vertical distance from the top
of the text to the current location).

Norm. Measures |  | $Elementary$ | $Advanced$ |
---|---|---|---|---
($/length$) | $\times 10^{n}$ | $\bar{X}$ | $\sigma$ | $r$ | $\bar{X}$ | $\sigma$ | $r$ | $Sign.$
Read time ($s$) | $\times 10^{2}$ | $3.4$ | $3.8$ | $0.10$ | $2.1$ | $1.9$ |  $0.14$ $\newmoon$ | $p<10^{-18}$
Regression num | $\times 10^{-3}$ | $5.9$ | $9.3$ | $-0.09$ | $7.1$ | $8.0$ | $-0.03$ | $p<10^{-3}$
Min speed ($px/ms$) | $\times 10^{-5}$ | $4.1$ | $11$ |  $-0.17$ $\newmoon$ | $3.3$ | $8.7$ | $0.01$ | -
Max speed ($px/ms$) | $\times 10^{-3}$ | $2.4$ | $1.7$ |  $-0.23$ $\newmoon$ | $1.7$ | $1.3$ | $-0.12$ | $p<10^{-12}$
Avg speed ($px/ms$) | $\times 10^{-4}$ | $7.9$ | $6.0$ |  $-0.35$ $\newmoon$ | $6.1$ | $4.6$ |  $-0.26$ $\newmoon$ | $p<10^{-14}$
Min acc ($px/ms^{2}$) | $\times 10^{-6}$ | $-5.9$ | $8.5$ | $0.12$ | $-4.4$ | $5.8$ | $0.08$ | $p<10^{-5}$
Max acc ($px/ms^{2}$) | $\times 10^{-6}$ | $-4.7$ | $8.7$ | $-0.11$ | $-3.4$ | $5.8$ | $-0.08$ | $p<10^{-3}$
Avg acc ($px/ms^{2}$) | $\times 10^{-7}$ | $-3.9$ | $13$ | $0.07$ | $-2.0$ | $6.7$ | $0.07$ | $p<10^{-3}$
$\newmoon$: $p<0.01$ $\bar{X}$: Mean value  $\sigma$: Standard deviation
Table 1: Interaction measures for $518$ participants across elementary and
advanced texts. Measures have been normalised according to text lengths. The
correlation ($r$) of measures with comprehension scores is presented (p-values
have been Bonferroni-corrected). The statistical significance of group
differences (Sign.) is calculated using a mixed-effects model.

###  3.1 Preprocessing

We collected responses from $600$ participants in total. However, the dataset
was initially preprocessed to mark entries where the participant had not
sufficiently engaged with reading the text. Given the size of the screen and
text lengths, scrolling was necessary in order to read the text. We only
included participants in our analysis if their scroll pattern indicated that
they had attempted to read the article. If no scrolling had been logged, or
the participant had not reached the half way point, the entry was not included
in the analysis. Removing these entries resulted in $518$ participants and
$1036$ articles marked as read.

### Interaction Measures

As a preliminary analysis, a range of interaction metrics were extracted from
the scrolling behaviour of participants. The computation of these follow the
standard of those used in prior scroll research Hinckley et al. (2002). These
measures are outlined below.

Elapsed time: The total reading time is recorded to produce a read time in
seconds.

Speed: The scroll speed ($s$) is calculated for each scroll interaction using
$s=d/t$, where $d$ represents the distance in pixels and $t$ the time taken in
$ms$. The average, minimum and maximum scroll speed are calculated.

Acceleration: Scrolling acceleration ($a$) is computed for each interaction.
This is calculated with the following formula: $a=(v-u)/t$ where $v$
represents the final scroll speed, $u$ the initial speed and $t$ denotes the
time taken in $ms$. The average, minimum and maximum scroll acceleration are
calculated.

Text regressions: Scrolling typically takes place in a linear vertical
fashion. Whilst reading, areas of text may require re-covering resulting in
upward scrolling actions. We count the number of times the participant
scrolled upwards to recover areas of text.

##  4 Results

### Readability Measures

Table 1 shows the mean values ($\bar{X}$) and the standard deviation
($\sigma$) of reading interactions normalised by text length. We additionally
present the correlation ($r$) of these with the comprehension scores of
participants. P-values have been Bonferroni-corrected to account for multiple
test conditions. The results show that three measures correlate significantly
with participant scores on elementary texts and two for advanced. For the
elementary articles, all speed measures correlate significantly with
participant scores. The negative correlation shows that the slower the speeds,
the higher the subsequent score. For advanced texts, the average scroll speed
also negatively correlates with score. Whereas the time taken to read the
article positively correlates with the subsequent score.

When considering the mean values across interactions, we see that all speed
and acceleration measures are lower for advanced texts. This finding shows
that, on average, the speed and acceleration of scrolling is slower on texts
that are more difficult. Additionally, the number of regressions is larger for
advanced texts than for beginner. Therefore, participants were more likely to
recover areas of text when it is at a higher level. Finally, the standard
deviation is larger for all measures on the elementary texts. This implies
that there is more variance in reading interaction styles when the text is at
a lower level.

We consider whether reading interactions differ significantly depending on the
level of the text. We calculate significance using Satterthwaite’s method
Kuznetsova et al. (2017), applied to a mixed-effects model that treats
participants and texts as crossed random effects.444Using R formula notation,
the model is: $measure\sim
readability+(readability|participant)+(readability|text)$. Tests were
performed using the lme4 and lmerTest R packages by Bates et al. (2014). All
measures are found to be statistically significant apart from the minimum
speed. The most significant measure ($p<10^{-18}$), is the normalised time a
participant spent on the text. Counter-intuitively, the time taken, on
average, is longer for the elementary texts than for advanced. However, when
we consider proficiency groups, participants with lower proficiency levels
took longer on elementary texts compared to advanced. This suggests that for
readers at a low level the text is too difficult to engage with, resulting in
the skipping of content. All acceleration measures, the maximum and average
speed as well as the number of regressions differ significantly depending on
whether text is advanced or elementary. These findings support the case that
there are different reading interactions for texts depending on their
complexity.

### Predicting Readability

System | $N$ | Precision | Recall | F-Score
---|---|---|---|---
Scroll only | 12 | 0.80 | 0.78 | 0.77
Vajjala-2018 | 155 | 0.88 | 0.85 | 0.84
Vajjala-2018 \+ Scroll | 160 | 0.92 | 0.88 | 0.88
Baseline (length + LR) | 6 | 0.91 | 0.87 | 0.88
Baseline + Traditional | 15 | 0.92 | 0.90 | 0.89
Baseline + Discourse | 24 | 0.81 | 0.82 | 0.80
Baseline + Syntactic | 24 | 0.87 | 0.85 | 0.84
Baseline + Psycholinguistic | 12 | 0.87 | 0.88 | 0.87
Baseline + Scroll | 18 | 0.98 | 0.97 | 0.96
Table 2: Table showing the total number of features per system ($N$) as well
as the precision, recall and f-score of models trained using 10-fold cross-
validation on the OneStopQA dataset.

We perform experiments to investigate whether scroll-based reading
interactions can be used to predict the level of a given text. As the texts in
our dataset have been read by multiple participants, we produce features by
taking the average of the statistically significant scroll measures (as
displayed under Sign. in Table 1). Prior research has shown that simple
methods of combining group interactions, such as averaging or majority voting,
can be highly robust Genre et al. (2013); Clemen (1989); Ertekin et al.
(2012). We produce two sets of scroll features using the number of
regressions, the max and average scroll speed and the max, min and average
scroll acceleration. One set is normalised by text length (Normalised) and the
other not ($\lnot$ Normalised). We then train a support vector machine (SVM)
to predict whether a text is advanced or elementary using these features. This
classifier was chosen as it has been shown to consistently yield better
results compared to other statistical models when predicting text readability
Kate et al. (2010). All reported results are obtained using stratified 10-fold
cross-validation.

Table 3 presents the f-score of the resulting models. To the best of our
knowledge, this is the first work aiming to predict readability using scroll-
based features. The best result ($0.77$) is achieved when using both sets of
scroll features. The $\lnot$ Normalised feature set performs better than
Normalised. This is likely due to the fact that these features contain signal
on the length of texts.

Scroll features | F-Score
---|---
All ($\lnot$Norm + Norm) | 0.77
$\lnot$ Normalised | 0.64
$\lnot$ Normalised | 0.61
Table 3: Readability prediction using aggregate interaction features.

Interestingly, filtering interactions by specific sub-groups can produce
better scores. The best results are achieved, using both sets of scroll
measures, when we only include the interactions of participants aged 25-34 –
resulting in an f-score of $0.81$. There are known differences in computer
interaction styles based on age Schneider et al. (2008). It therefore follows
that aggregate measures from an audience with a more consistent interaction
style would result in better features.

The use of scroll features alone produces an f-score of $0.77$. We compare our
result to a state-of-the-art readability system by Vajjala and Lučić (2018)
referred to as Vajjala-2018. This system uses a multilayer perceptron
classifier and has been shown to outperform BERT-based approaches on the
OneStopEnglish dataset Martinc et al. (2019). The system relies on 155 hand-
crafted features which are grouped into six categories: traditional metrics,
word features, psycholinguistic, lexical richness, syntactic and discourse
features. The Vajjala-2018 system outperforms the use of scroll features alone
when trained to predict if texts are advanced or elementary on our dataset.
However, we are able to improve the f-score of this state-of-the-art system by
$4\%$ when adding aggregate scroll features. This is an interesting result as
it shows that 1) these features are highly complementary to existing
readability systems and 2) these features represent an aspect of textual
complexity not covered by the current $155$ features. This may be due to
scroll-based features capturing a notion of conceptual complexity which plays
a vital role in text understanding and maintaining a reader’s interest Štajner
and Hulpus (2020).

We also investigate how scroll-interaction features compare to classical
readability feature sets (an overview of feature sets is included in Appendix
A). To do this, we initially create a highly competitive baseline by training
an SVM using the length of text and measures of lexical richness. Due to the
nature of the OneStopEnglish dataset, length is an extremely informative
feature. This is because, on average, the word length for advanced texts
($915$) is almost always higher than for the simplified elementary texts
($599$). Vocabulary knowledge, including lexical diversity and richness, are
principal components in reading comprehension Collins-Thompson (2014). The
importance of lexical richness has been investigated from the perspective of
second language acquisition Lu (2012), as well as in readability systems
(Vajjala and Meurers, 2012; Vajjala and Lučić, 2018; Xia, 2019). We opted to
use features in our baseline system that were highly informative but would not
require extensive text processing. A key advantage of lexical richness
measures is that they are a function of the number of word types (T) and the
total text length (N). Therefore, they can be calculated quickly and are
robust to noisy or broken text. The specific measures of lexical richness we
include are type-token ratio (TTR), which is the ratio of the number of unique
word tokens to the total number of word tokens in a text, Root TTR
$(T/\sqrt{N})$, Corrected TTR $(T/\sqrt{2N})$, Bilogarithmic TTR ($LogT/LogN$)
and Uber Index ($Log^{2}T/Log(N/T)$). The resulting baseline approach
outperforms the Vajjala-2018 system. However, it should be noted that the
Vajjala-2018 system was originally designed to differentiate between
elementary, intermediate and advanced texts which is a more nuanced and
difficult task.

By performing ablation tests, we see how scroll-interaction features compare
to classical sets of readability features when combined with a strong baseline
system (Length + LR). We find that adding sets of discourse, syntactic or
psycholinguistic features degrades classifier performance. The addition of
traditional measures, such as Flesch-Kincaid score, produces a small
improvement. However, adding aggregate scroll features produces a marked
improvement on the baseline resulting in an f-score of $0.96$. The addition of
scroll-interaction features improves performance beyond all other classical
feature combinations. A key advantage to this approach is that the result is
achieved with only $18$ features which are highly robust to noisy text.

### Subjective Readability

Readability approaches typically produce an objective numerical score which
often corresponds to a suggested level Fry (2002). In our previous
experiments, we predict the readability of text as defined by such preordained
levels (i.e., elementary or advanced). However, the readability of text can
also be defined in a more subjective and idiosyncratic manner. Such techniques
have been referred to as levelling, and are similar to readability in that
they determine difficulty but are subjective Clay (1991). Levelling integrates
a reader’s background and experience with objective readability by
understanding what contributes to readability for differing audiences. In this
section, we investigate what reading-interactions could tell us about an
individual’s text comprehension or L1.

In our study, we record comprehension scores to evaluate the understanding and
readability of text for individuals. Participants are asked three questions
per article assessing the reader’s understanding of what they have read.
Figure 4 shows the average score of participants grouped by their self-
reported proficiency. We see that the average score decreases in line with the
reported proficiency.

In the fashion of levelling, we investigate the correlation of features with
the comprehension scores of audiences. We observe whether the feature
correlations, at both the text and interaction level, vary depending on the
first language (L1) of the group. We focus on English (n = $350$) and Tamil (n
= $101$). Tamil was chosen as it was the second most represented L1 after
English.

We find there is a statistically significant association between the average
scroll speed and scores of participants for both Tamil and English
($p<10^{-4}$). The correlation is negative, indicating that the faster the
average speed of reading, the lower the subsequent score. When considering the
English L1 audience, the only statistically significant correlation with score
was the average reading speed. However, for Tamil readers there were three
statistically significant correlates with ($p<10^{-4}$), in order of strength
these were: the mean age of acquisition for vocabulary in the text (AoA) as
reported by Kuperman et al. (2012); the average reading speed and the mean AoA
of word lemmas. The negative correlation shows that the higher the AoA of
vocabulary, the lower the score for this group. This finding supports prior
work emphasising the importance of AoA as a factor when simplifying texts for
second language learners Crossley et al. (2007).

Figure 4: Bar chart showing the average scores (out of three) for texts across
proficiency levels.

Figure 5: Histograms showing the density of average scroll speed for Tamil and
English as well as across participant scores.

The average scrolling speed correlates significantly with scores for both
English and Tamil L1 audiences. We investigate whether the distribution of
average reading speeds differs based on the reader group. Figure 5 presents a
histogram of the average reading speeds. In both groups, we observe a
positively right skewed distribution, as speed is a vector quantity and a
range boundary occurs at 0. When comparing speed distributions, the mean and
median speeds are higher for English ($0.50$, $0.47$) than Tamil ($0.46$,
$0.37$). The higher median and mean measures for English interactions are
likely due to a higher familiarity yielding a faster reading speed. The L1
English group additionally achieve a higher score on average ($1.18$) than the
Tamil audience ($0.98$), despite the faster scroll speed.

In Figure 5, we consider how the distribution of average reading speeds varies
for participants according to the score they attained. The higher the score,
on average, the lower the scroll (reading) speed. For the lower scores, we see
that the distribution has a wider spread. According to the self-reported
proficiency ratings, the advanced texts would likely have been too difficult
for a proportion of participants. A higher scrolling speed indicates that a
reader is skipping content without properly reading, perhaps due to the level
being too high for the reader to sufficiently engage with. Finally, there is a
statistically significant positive correlation between the average reading
speed of a participant and their reported proficiency ($p<10^{-4}$), further
supporting the notion that reading interactions vary based on ability.

Being able to understand how on-device reading interactions vary according to
an individuals’ L1 and comprehension is incredibly useful. Such information
could be used in text simplification systems, or in a ‘levelling’ manner to
match the appropriate level of text to a given reader. Such applications are
especially useful for individuals who are learning a language.

##  5 Conclusions

To conclude, the use of scroll features for judging readability has numerous
benefits. Such measures are language agnostic, unobtrusive and are robust to
noisy text. Furthermore, implicit user feedback allows an insight into
readability at an individual level, thereby allowing for a more inclusive and
personalisable assessment. We present a $518$ participant study to investigate
the impact of text readability on reading interactions. In this paper, we
confirm that there are statistically significant differences in the way that
readers interact with advanced and elementary texts, and that the
comprehension scores of individuals correlate with specific measures of
scrolling interaction. We demonstrate that, even with a simple model,
aggregate scroll interactions can be used to predict readability. Moreover, we
show that individual scroll behaviour can provide an insight into the
subjective readability for an individual. Finally, we improve a state-of-the-
art readability classifier with the integration of scroll-interaction
features, demonstrating that interaction features are highly complementary to
traditional linguistic approaches. In future work, we will focus on
investigating which aspects of readability scroll-based measures index.

## References

  * Acheson et al. (2008) Daniel J Acheson, Justine B Wells, and Maryellen C MacDonald. 2008.  New and updated tests of print exposure and reading abilities in college students.  _Behavior research methods_ , 40(1):278–289.
  * Agić et al. (2016) Željko Agić, Anders Johannsen, Barbara Plank, Héctor Martínez Alonso, Natalie Schluter, and Anders Søgaard. 2016\.  Multilingual projection for parsing truly low-resource languages.  _Transactions of the Association for Computational Linguistics_ , 4:301–312.
  * Aluisio et al. (2010) Sandra Aluisio, Lucia Specia, Caroline Gasperin, and Carolina Scarton. 2010.  Readability assessment for text simplification.  In _Proceedings of the NAACL HLT 2010 Fifth Workshop on Innovative Use of NLP for Building Educational Applications_ , pages 1–9.
  * Bates et al. (2014) Douglas Bates, Martin Mächler, Ben Bolker, and Steve Walker. 2014.  Fitting linear mixed-effects models using lme4.  _arXiv preprint arXiv:1406.5823_.
  * Berzak et al. (2020) Yevgeni Berzak, Jonathan Malmaud, and Roger Levy. 2020.  Starc: Structured annotations for reading comprehension.  _arXiv preprint arXiv:2004.14797_.
  * Chall (1958) Jeanne Sternlicht Chall. 1958.  _Readability: An appraisal of research and application_.  34\. Ohio State University.
  * Chen et al. (2021) Robert Chen, Roger Levy, and Tiwalayo Eisape. 2021.  On factors influencing typing time: Insights from a viral online typing game.  In _Proceedings of the Annual Meeting of the Cognitive Science Society_ , volume 43.
  * Clay (1991) Marie M Clay. 1991.  _Becoming literate: The construction of inner control_.  Heinemann Educational Books.
  * Claypool et al. (2001) Mark Claypool, Phong Le, Makoto Wased, and David Brown. 2001.  Implicit interest indicators.  In _Proceedings of the 6th international conference on Intelligent user interfaces_ , pages 33–40.
  * Clemen (1989) Robert T Clemen. 1989.  Combining forecasts: A review and annotated bibliography.  _International journal of forecasting_ , 5(4):559–583.
  * Collins-Thompson (2013) Kevyn Collins-Thompson. 2013.  Enriching the web by modeling reading difficulty.  In _ESAIR_ , pages 3–4.
  * Collins-Thompson (2014) Kevyn Collins-Thompson. 2014.  Computational assessment of text readability: A survey of current and future research.  _ITL-International Journal of Applied Linguistics_ , 165(2):97–135.
  * Collins-Thompson and Callan (2004) Kevyn Collins-Thompson and James P Callan. 2004.  A language modeling approach to predicting reading difficulty.  In _Proceedings of the Human Language Technology Conference of the North American Chapter of the Association for Computational Linguistics: HLT-NAACL 2004_ , pages 193–200.
  * Crossley et al. (2007) Scott A Crossley, Philip M McCarthy, and Danielle S McNamara. 2007.  Discriminating between second language learning text-types.  In _FLAIRS Conference_ , pages 205–210.
  * Crossley et al. (2017) Scott A Crossley, Stephen Skalicky, Mihai Dascalu, Danielle S McNamara, and Kristopher Kyle. 2017.  Predicting text comprehension, processing, and familiarity in adult readers: New approaches to readability formulas.  _Discourse Processes_ , 54(5-6):340–359.
  * Dale and Chall (1949) Edgar Dale and Jeanne S Chall. 1949.  The concept of readability.  _Elementary English_ , 26(1):19–26.
  * Ertekin et al. (2012) Seyda Ertekin, Haym Hirsh, and Cynthia Rudin. 2012.  Learning to predict the wisdom of crowds.  _arXiv preprint arXiv:1204.3611_.
  * Feng et al. (2009) Lijun Feng, Noémie Elhadad, and Matt Huenerfauth. 2009.  Cognitively motivated features for readability assessment.  In _Proceedings of the 12th Conference of the European Chapter of the ACL (EACL 2009)_ , pages 229–237.
  * Feng et al. (2010) Lijun Feng, Martin Jansche, Matt Huenerfauth, and Noémie Elhadad. 2010.  A comparison of features for automatic readability assessment.
  * François and Miltsakaki (2012) Thomas François and Eleni Miltsakaki. 2012.  Do nlp and machine learning improve traditional readability formulas?  In _Proceedings of the First Workshop on Predicting and Improving Text Readability for target reader populations_ , pages 49–57.
  * Fry (2002) Edward Fry. 2002.  Readability versus leveling.  _The Reading Teacher_ , 56(3):286–291.
  * Genre et al. (2013) Véronique Genre, Geoff Kenny, Aidan Meyler, and Allan Timmermann. 2013.  Combining expert forecasts: Can anything beat the simple average?  _International Journal of Forecasting_ , 29(1):108–121.
  * Golovchinsky et al. (1999) Gene Golovchinsky, Morgan N Price, and Bill N Schilit. 1999.  From reading to retrieval: freeform ink annotations as queries.  In _Proceedings of the 22nd annual international ACM SIGIR conference on Research and development in information retrieval_ , pages 19–25.
  * González-Garduno and Søgaard (2018) Ana V González-Garduno and Anders Søgaard. 2018.  Learning to predict readability using eye-movement data from natives and learners.  In _Thirty-Second AAAI Conference on Artificial Intelligence_.
  * Graesser et al. (2011) Arthur C Graesser, Danielle S McNamara, and Jonna M Kulikowich. 2011.  Coh-metrix: Providing multilevel analyses of text characteristics.  _Educational researcher_ , 40(5):223–234.
  * Heilman et al. (2007) Michael Heilman, Kevyn Collins-Thompson, Jamie Callan, and Maxine Eskenazi. 2007\.  Combining lexical and grammatical features to improve readability measures for first and second language texts.  In _Human Language Technologies 2007: The Conference of the North American Chapter of the Association for Computational Linguistics; Proceedings of the Main Conference_ , pages 460–467.
  * Heilman et al. (2008) Michael Heilman, Kevyn Collins-Thompson, and Maxine Eskenazi. 2008.  An analysis of statistical models and features for reading difficulty prediction.  In _Proceedings of the third workshop on innovative use of NLP for building educational applications_ , pages 71–79.
  * Hinckley et al. (2002) Ken Hinckley, Edward Cutrell, Steve Bathiche, and Tim Muss. 2002.  Quantitative analysis of scrolling techniques.  In _Proceedings of the SIGCHI conference on Human factors in computing systems_ , pages 65–72.
  * Kate et al. (2010) Rohit Kate, Xiaoqiang Luo, Siddharth Patwardhan, Martin Franz, Radu Florian, Raymond Mooney, Salim Roukos, and Chris Welty. 2010.  Learning to predict readability using diverse linguistic features.  In _Proceedings of the 23rd International Conference on Computational Linguistics (Coling 2010)_ , pages 546–554.
  * Kelly and Teevan (2003) Diane Kelly and Jaime Teevan. 2003.  Implicit feedback for inferring user preference: a bibliography.  In _Acm Sigir Forum_ , volume 37, pages 18–28. ACM New York, NY, USA.
  * Kincaid et al. (1975) J Peter Kincaid, Robert P Fishburne Jr, Richard L Rogers, and Brad S Chissom. 1975\.  Derivation of new readability formulas (automated readability index, fog count and flesch reading ease formula) for navy enlisted personnel.  Technical report, Naval Technical Training Command Millington TN Research Branch.
  * Klare et al. (1963) George Roger Klare et al. 1963.  Measurement of readability.
  * Kuperman et al. (2012) Victor Kuperman, Hans Stadthagen-Gonzalez, and Marc Brysbaert. 2012.  Age-of-acquisition ratings for 30,000 english words.  _Behavior research methods_ , 44(4):978–990.
  * Kuznetsova et al. (2017) Alexandra Kuznetsova, Per B Brockhoff, Rune HB Christensen, et al. 2017.  lmertest package: tests in linear mixed effects models.  _Journal of statistical software_ , 82(13):1–26.
  * Li et al. (2019) Qisheng Li, Meredith Ringel Morris, Adam Fourney, Kevin Larson, and Katharina Reinecke. 2019.  The impact of web browser reader views on reading speed and user experience.  CHI ’19, page 1–12, New York, NY, USA. Association for Computing Machinery.
  * Lu (2012) Xiaofei Lu. 2012.  The relationship of lexical richness to the quality of esl learners’ oral narratives.  _The Modern Language Journal_ , 96(2):190–208.
  * Martinc et al. (2019) Matej Martinc, Senja Pollak, and Marko Robnik-Šikonja. 2019.  Supervised and unsupervised neural approaches to text readability.  _arXiv preprint arXiv:1907.11779_.
  * Morris et al. (2018) Meredith Ringel Morris, Adam Fourney, Abdullah Ali, and Laura Vonessen. 2018.  Understanding the needs of searchers with dyslexia.  CHI ’18, page 1–12, New York, NY, USA. Association for Computing Machinery.
  * Pitler and Nenkova (2008) Emily Pitler and Ani Nenkova. 2008.  Revisiting readability: A unified framework for predicting text quality.  In _Proceedings of the 2008 conference on empirical methods in natural language processing_ , pages 186–195.
  * Schneider et al. (2008) Nicole Schneider, Janet Wilkes, Morten Grandt, and Christopher M Schlick. 2008.  Investigation of input devices for the age-differentiated design of human-computer interaction.  In _Proceedings of the Human Factors and Ergonomics Society Annual Meeting_ , volume 52, pages 144–148. SAGE Publications Sage CA: Los Angeles, CA.
  * Si and Callan (2001) Luo Si and Jamie Callan. 2001.  A statistical model for scientific readability.  In _Proceedings of the tenth international conference on Information and knowledge management_ , pages 574–576.
  * Štajner and Hulpus (2020) Sanja Štajner and Ioana Hulpus. 2020.  When shallow is good enough: Automatic assessment of conceptual text complexity using shallow semantic features.  In _Proceedings of The 12th Language Resources and Evaluation Conference_ , pages 1414–1422.
  * Vajjala and Lučić (2018) Sowmya Vajjala and Ivana Lučić. 2018.  Onestopenglish corpus: A new corpus for automatic readability assessment and text simplification.  In _Proceedings of the thirteenth workshop on innovative use of NLP for building educational applications_ , pages 297–304.
  * Vajjala and Meurers (2012) Sowmya Vajjala and Detmar Meurers. 2012.  On improving the accuracy of readability classification using insights from second language acquisition.  In _Proceedings of the seventh workshop on building educational applications using NLP_ , pages 163–173.
  * White et al. (2002) Ryen W White, Ian Ruthven, and Joemon M Jose. 2002.  Finding relevant documents using top ranking sentences: an evaluation of two alternative schemes.  In _Proceedings of the 25th annual international ACM SIGIR conference on Research and development in information retrieval_ , pages 57–64.
  * Xia (2019) Menglin Xia. 2019.  _Text readability and summarisation for non-native reading comprehension_.  Ph.D. thesis, University of Cambridge.
  * Xia et al. (2016) Menglin Xia, Ekaterina Kochmar, and Ted Briscoe. 2016.  Text readability assessment for second language learners.  In _Proceedings of the 11th Workshop on Innovative Use of NLP for Building Educational Applications_ , pages 12–22.
  * Zakaluk and Samuels (1988) Beverly L Zakaluk and S Jay Samuels. 1988.  _Readability: Its Past, Present, and Future._ ERIC.
  * Zou et al. (2019) Yixin Zou, Shawn Danino, Kaiwen Sun, and Florian Schaub. 2019.  You ‘might’ be affected: An empirical analysis of readability and usability issues in data breach notifications.  In _Proceedings of the 2019 CHI Conference on Human Factors in Computing Systems_ , CHI ’19, page 1–14, New York, NY, USA. Association for Computing Machinery.

##  6 Appendix

### A - Traditional Readability Feature Sets

The readability feature sets for the OneStopEnglish dataset were provided by
(Vajjala and Lučić., 2018)555https://zenodo.org/record/1219041. Please refer
to their work for further details. An overview of features is provided:

#### Traditional

  * •

Avg. number of characters per word

  * •

Avg. number of syllables per word

  * •

Avg. sentence length

  * •

Flesch-Kincaid score

  * •

Coleman-Liau readability formula

  * •

SMOG grade

#### Discourse

  * •

Word overlap features - content word, noun, stem and argument overlap at local
(between adjacent sentences) and global (between any two sentences in a text)
levels.

  * •

Entity transitions - features based on the transitions between the syntactic
roles of entities.

  * •

Co-reference chains - features based on noun phrases, nouns and pronouns and
determiner usage.

  * •

Num. of referential expressions.

  * •

Num. of discourse and non-discourse connectives and all connectives

#### Syntactic

  * •

Num. NPs per sentence (NumNP)

  * •

Num. VPs per sentence (NumVP)

  * •

Num. PPs per sentence (NumPP))

  * •

Avg. length of a NP (NPSize)

  * •

Num. Dependent Clauses per sentence (NumDC)

  * •

Num. Complex-T units per sentence (NumCT)

  * •

Num. Co-ordinate Phrases per sentence (CoOrd)

  * •

Num. SBARs per sentence (NumSBAR)

  * •

Avg. Parse Tree Height (TreeHeight)

  * •

Avg. length of a VP (VPSize)

  * •

Avg. length of a PP (PPSize)

  * •

Mean length of clause (MLC)

  * •

Mean length of a sentence (MLS)

  * •

Mean length of T-unit (MLT)

  * •

Num. of Clauses per Sentence (C/S)

  * •

Num. of T-Units per sentence (T/S)

  * •

Num. of Clauses per T-unit (C/T)

  * •

Num. of Complex-T-Units per T-unit (CT/T)

  * •

Dependent Clause to Clause Ratio (DC/C)

  * •

Dependent Clause to T-unit Ratio (DC/T)

  * •

Co-ordinate Phrases per Clause (CP/C)

  * •

Co-ordinate Phrases per T-unit (CP/T)

  * •

Complex Nominals per Clause (CN/C) – Complex Nominals per T-unit (CN/T)

  * •

Verb phrases per T-unit (VP/T)

#### Psycholinguistic

Norms from MRC which were compiled by Gilhooly and Logie (1980) for 1903
unique words including:

  * •

Avg. word age of acquisition

  * •

Avg. word Familiarity

  * •

Avg. word concreteness

  * •

Avg. word imagability

  * •

Avg. word meaningfulness


# Radiative thermal switch exploiting hyperbolic surface phonon polaritons

Annika Ott1, Yang Hu2,3,4 , Xiao-Hu Wu4, Svend-Age Biehs1 s.age.biehs@uni-
oldenburg.de 1 Institut für Physik, Carl von Ossietzky Universität, D-26111
Oldenburg, Germany 2 Basic Research Center, School of Power and Energy,
Northwestern Polytechnical University, Xi’an, 710064, P.R. China 3 Center of
Computational Physics and Energy Science, Yangtze River Delta Research
Institute of NPU, Northwestern Polytechnical University, Taicang, Jiangsu,
215400, P.R. China 4 Shandong Institute of Advanced Technology, Jinan 250100,
P.R. China

(August 27, 2024)

###### Abstract

We study the radiative heat flux between two nanoparticles in close vicinity
to the natural hyperbolic material hBN with its optical axis oriented parallel
to the interface. We show that the heat flux between the nanoparticles can be
efficiently modulated when rotating the nanoparticles due to the coupling to
the highly directional hyperbolic surface modes in hBN. Finally, we discuss
the thickness and distance dependence of this effect.

##  I Introduction

Recently, several possibilities for the manipulation of near-field radiative
heat transport due to many-body interactions have been highlighted: It could
be shown that the radiative heat flux (HF) between two objects can be
tremendously enhanced by the interaction with an environment which has
interfaces supporting long-range surface modes Saaskilathi2014 ; Asheichyk2017
; DongEtAl2018 ; paper_2sic ; HeEtAl2019 . Such an enhancement can also be
observed for the HF through plasmonic or hyperbolic films MessinaEtAl2012 ;
MessinaEtAl2016 ; ZhangEtAl2019b ; McSherryLenert2020 due to the transmission
of large wave-vector modes through such structures. Recently, non-reciprocal
media such as magneto-optical materials or Weyl metals have been studied in
many-body assemblies because they exhibit, on the one hand, interesting
effects like a persistent heat current zhufan ; zhufan2 , giant magneto-
resistance Latella2017 ; Cuevas ; HeEtAl2020 ; Cuevas2 ; Song , a Hall and
anomalous Hall effect for thermal radiation hall ; Ahall as well as a circular
near-field HF which is closely connected spin and angular momenta of thermal
radiation Ott2018 ; Zubin2019 ; Khandekar . On the other hand, such non-
reciprocal media allow for actively controlling the directionality of surface
modes and therefore can be used to control the strength of HF between two
objects via the non-reciprocal surface modes of the environment which can
result in a strong HF rectification Ott2019 ; Ott2020 . Another possibility is
to use the intrinsic anisotropy of the environment as encountered in uni-axial
materials as shown for the HF between two nanoparticles (NPs) in the vicinity
of a 2D meta-surface made of graphene strips ZhangEtAl2019 ; graphene also
offers the possibility for active thermal switching by electrical biasing
YuEtAl2017 ; IlicEtAl2018 which is much more efficient than electrical
switching with ferroelectric materials HuangEtAl2014 . A review of these
effects and other recent developments in the field of many-body effects for
near-field thermal radiation can be found in Ref. BiehsEtAl2021 .

In this work we propose a thermal switching effect for the HF between two hBN
NPs in close vicinity to a hBN substrate. hBN is a natural hyperbolic material
Narimanov which supports like black phosphorus or periodically ordered
graphene or graphite structures directional hyperbolic surface modes
LiuEtAl2019 ; ShenEtAl2018 ; LiuXuan2016 ; Wu2021 ; Wu2021c . Due to its
intrinsic anisotropy the radiative HF between two hBN sheets can be modulated
by the relative twisting of their optical axes LiuEtAl2017 ; Wu2021b an effect
which is similar to the HF modulation between grating structures BiehsEtAl2011
. In our work, we will show that the strong directionality of the energy or
heat flow of the hyperbolic surface phonon polaritons (HSPP) in an hBN
substrate can be used to modulate the radiative HF between two NPs by a huge
amount when rotating the substrate or the NPs. We find that the modulation
contrast for the full HF is about almost 1500 which is in stark contrast to
the modulation contrast found between two hBN films which is about 5.36 at its
maximum LiuEtAl2017 ; Wu2021b and between two gratings where it is about 10 at
its maximum BiehsEtAl2011 .

Figure 1: Sketch of the considered configuration together with a visualization
of the two HF channels: direct radiative heat exchange and heat exchange via
modes of the substrate like wave-guide modes and surface modes. The angle
$\phi$ is the rotation angle between the optical axis and the axis connecting
the two NPs.

##  II Setup and formalism

We consider a setup as sketched in Fig. 1 where two spherical hBN NPs (labeled
with 1 and 2) are at a distance $z$ from a planar hBN substrate. The two NPs
have temperatures $T_{1}$ and $T_{2}$, respectively, and the inter-particle
distance is $d$. When omitting the thermal emission to the background
Domingues ; Chapuis ; Ott2020 , the net reveiced/emitted power for particle
$i=1,2$ can be written as ($j=1,2;j\neq i$)

$P_{i}=\int_{0}^{\infty}\\!\\!\frac{{\rm
d}\omega}{2\pi}\,P_{i,\omega}=3\int_{0}^{\infty}\\!\\!\frac{{\rm
d}\omega}{2\pi}\,\bigl{[}\Theta(T_{j})-\Theta(T_{i})\bigr{]}\mathcal{T}_{ji}$
(1)

where $\Theta(T)=\hbar\omega/(\exp(\hbar\omega/k_{B}T)-1)$ is the Bose-
Einstein function; $\hbar$ is the reduced Planck constant and $k_{B}$ is the
Boltzmann constant. $\mathcal{T}_{ji}$ is the transmission coefficients for
the power exchanged between the NPs and is given by EkerothEtAl2017

$\displaystyle\mathcal{T}_{ji}$ $\displaystyle=\frac{4}{3}k_{0}^{4}{\rm
Tr}\bigl{[}\uuline{\chi}\mathds{G}_{ij}\uuline{\chi}\mathds{G}_{ij}^{\dagger}\bigr{]},$
(2)

where we have introduced the wave number in vacuum $k_{0}=\omega/c$. When
neglecting the radiative correction we can write the response function of the
NPs $\uuline{\chi}$ in terms of the polarizability tensor $\uuline{\alpha}$ as
EkerothEtAl2017

$\displaystyle\uuline{\chi}=\frac{1}{2{\rm
i}}\bigl{(}\uuline{\alpha}-\uuline{\alpha}^{\dagger}\bigr{)}.$ (3)

Furthermore, $\mathds{G}_{ij}=\mathds{G}(\mathbf{r}_{i},\mathbf{r}_{j})$ are
the Green functions at positions $\mathbf{r}_{i/j}$ of particle $i$ or $j$.
The detailed definition can be found in App. B. Finally, for spherical NPs
having a radius $R$ the polarizability in dipole approximation (quasi-static
limit) is given by LakhtakiaEtAl1991

$\underline{\underline{\alpha}}=4\pi
R^{3}(\underline{\underline{\epsilon}}-\mathds{1})(\underline{\underline{\epsilon}}+2\mathds{1})^{-1}.$
(4)

As a consequence, the heat exchange scales like $R^{6}$ in the distance regime
$d,z>4R$ where the dipole model can be applied.

##  III HVPP and HSPP in hBN

Figure 2: Regions in the $q_{x}$-$q_{y}$ plane where HVPP and HSPP can exist
for a) $\omega=1.5\times 10^{14}$rad/s (type I) and b) $\omega=2.9\times
10^{14}$rad/s (type II) when the optical axis (OA) is along the y direction.
c) and d) show the corresponding ${\rm Im}(r_{pp})$. The green dashed lines
are the iso-frequency lines from Eq. (8) dividing the $q_{x}$-$q_{y}$ plane
into the regions with HVPP and HSPP. The magenta dashed lines are the HVPP and
HSPP solutions of Eq. (10) and the white dashed lines are the solutions of Eq.
(9).

In the following we assume that the NPs and the substrate are made of hBN and
therefore are natural hyperbolic materials Narimanov . Furthermore, we assume
that the optical axis is along the y direction so that the permittivity tensor
can be written as

$\underline{\underline{\rm\epsilon}}={\rm
diag}(\epsilon_{\perp},\epsilon_{\parallel},\epsilon_{\perp})$ (5)

with the permittivity along/perpendicular to the optical axis LiuXuan2016

$\epsilon_{\perp,\parallel}=\epsilon_{\infty\perp,\parallel}\biggl{(}1+\frac{\omega_{L\perp,\parallel}^{2}-\omega_{T\perp,\parallel}^{2}}{\omega_{T\perp,\parallel}^{2}-\omega_{\perp,\parallel}^{2}-{\rm
i}\gamma_{\perp,\parallel}\omega}\biggr{)}$ (6)

with $\epsilon_{\perp}=2.95$, $\omega_{L\perp}=1610$cm-1 ,
$\omega_{T\perp}=1610$cm-1, $\gamma_{\perp}=5$cm-1 and
$\epsilon_{\perp}=4.87$, $\omega_{L\perp}=830$ cm-1, $\omega_{T\perp}=780$
cm-1 and $\gamma_{\perp}=4$cm-1. As already verified in several other works
LiuXuan2016 ; LiuEtAl2017 ; Wu2021 ; Wu2021b ; Wu2021c there is a hyperbolic
band of type I (${\rm Re}(\epsilon_{\parallel})<0,{\rm
Re}(\epsilon_{\perp})>0$) in the frequency range $1.47-1.56\times 10^{14}$
rad/s and a hyperbolic band of type II (${\rm Re}(\epsilon_{\parallel})>0,{\rm
Re}(\epsilon_{\perp})<0$) in the frequency range $2.58-3.03\times 10^{14}$
rad/s.

Due to its natural anisotropy hBN can support hyperbolic volume phonon
polaritons (HVPP) as well as HSPP for the extra-ordinary modes inside the uni-
axial material. The HVPP are the propagating wave solutions for waves
traveling inside hBN within the hyperbolic frequency bands. When considering
thin films of hBN these waves can form Fabry-Pérot-like standing waves within
the thin film. On the other hand, the HSPP are evanescent wave solutions of
the wave equation within the hyperbolic frequency bands which can travel
within the interface but have exponentially decaying electromagnetic fields
perpendicular to the interface. Therefore the HSPP are confined to the
interfaces of hBN and vacuum. For thin films of hBN the HSPP of the two
interfaces can couple leading to hybridized HSPPs. The $k_{x}$-$k_{y}$ range
where either HVPP or HSPP can exist, can therefore be distinguished by the $z$
component of the wave vector for the considered extra-ordinary modes(eom)
LiuEtAl2017 ; Wu2021 , which determines whether the fields in z-direction are
exponentially damped or not. In our case it is given by LiuEtAl2017 ;
LiEtAl2017 ; WuEtAl2019

$k_{z,\rm
eom}=\sqrt{k_{0}^{2}\epsilon_{\parallel}-k_{x}^{2}-k_{y}^{2}\epsilon_{\parallel}/\epsilon_{\perp}}.$
(7)

When this component is purely imaginary for real-valued permittivities, then
the corresponding modes would be evanescent within the uni-axial materials
which are the HSPP if $k_{x},k_{y}>k_{0}$. In the case that it is purely real-
valued for real-valued permittivities the corresponding modes are propagating
HVPP. Both regions are visualized in both hyperbolic bands in Fig. 2a) and b).
The HVPP only exist above and below the isofrequency lines which follow from
the condition $k_{z,\rm eom}=0$, i.e.

$\frac{k_{x}^{2}}{k_{0}^{2}\epsilon_{\parallel}}+\frac{k_{y}^{2}}{k_{0}^{2}\epsilon_{\perp}}=1,$
(8)

whereas the HSPP can only exist on the left and right of these iso-frequency
lines fulfilling the dispersion relation LiuEtAl2017

$\frac{k_{y}^{2}}{k_{0}^{2}}(1-\epsilon_{\parallel}\epsilon_{\perp})+\frac{k_{x}^{2}}{k_{0}^{2}}(1-\epsilon_{\perp}^{2})=-\epsilon_{\parallel}\epsilon_{\perp}.$
(9)

Note that the HSPP do not exist in hBN for the type I case Wu2021 . For thin
films similar but more elaborate dispersion relations can be derived
AlvarezEtAl2019 . In the case that $k_{x},k_{y}\gg k_{0}$ it reads
($\l=0,1,2,\ldots$)

$\sqrt{k_{x}^{2}+k_{y}^{2}}=\frac{2\rho}{t}\biggl{[}\arctan\biggl{(}\frac{\rho}{\epsilon_{\perp}}\biggr{)}+\frac{\pi}{2}\l\biggr{]}$
(10)

with $\rho={\rm
i}\sqrt{\epsilon_{\perp}(k_{x}^{2}+k_{y}^{2})/(\epsilon_{\perp}k_{x}^{2}+\epsilon_{\parallel}k_{y}^{2})}$.
Note, that this relation holds for HVPP and HSPP. The corresponding results
are plotted in Fig. 2c) and d).

Due to the fact that the optical axis is within the interface, there are
depolarization effects which means that an incoming s-polarized wave can be
scattered into a p-polarized wave and vice versa. Therefore the reflection
coefficients $r_{sp}$ and $r_{ps}$ are in general non-zero and therefore the
reflection tensor

$\mathds{R}=\begin{pmatrix}r_{ss}&r_{sp}\\\ r_{ps}&r_{pp}\end{pmatrix}$ (11)

is not diagonal anymore. However, the near-field coupling of the two NPs via
the HSPP of the substrate is mainly given by the component $r_{pp}$ of the
reflection tensor describing the reflection of an incoming p-polarized wave
into a p-polarized wave so that we can understand the underlying physics by
focusing on that component. In Fig. 2 we show ${\rm Im}(r_{pp})$ for an hBN
substrate of thickness $t=500$nm with optical axis in y direction for two
different frequencies from the hyperbolic type I and II bands within the
$q_{x}$-$q_{y}$ plane where $q_{x/y}=k_{x/y}/k_{0}$. Furthermore, we have
added the corresponding dispersion relations of the HVPP and HSPP. Note, that
the dispersion relation for semi-infinite materials coincides with the
corresponding dispersion relation for a thin film when $q_{x},q_{y}\gg 1$.

##  IV Numerical results

Figure 3: ${\rm Im}(\alpha_{x})$ and ${\rm Im}(\alpha_{y})$ for hbN NPs with
optical axis in y direction normalized to the particle volume $4\pi R^{3}/3$.
Note that ${\rm Im}(\alpha_{z})={\rm Im}(\alpha_{x})$. The horizontal lines
mark the edges of the two hyperbolic frequency bands.

Now, when bringing the two hBN NPs in close vicinity to the substrate they can
couple to the HVPP and HSPP and thus the heat flow between the two NPs can be
tremendously enhanced, in particular, via the surface modes Saaskilathi2014 ;
Asheichyk2017 ; DongEtAl2018 ; paper_2sic . In order to see which of the two
hyperbolic frequency bands is more relevant for the heat transport we show in
Fig. 3 the imaginary part of the components of the polarizability which is
proportional to the absorptivity of the NPs. It can be observed that the
component parallel to the optical axis has a resonance in the low-frequency
hyperbolic band of type I and the other two components have a resonance in the
high-frequency hyperbolic band of type II. This second resonance at
$\omega=2.906\times 10^{14}$rad/s is much stronger than the first one so that
the HF between the two NPs will be dominated by the type II band. Therefore,
mainly the HSPP in this high-frequency band are responsible for the radiative
coupling between the NPs.

Figure 4: Imaginary part of $r_{pp}$ at $\omega=2.9\times 10^{14}\,{\rm
rad/s}$ as in Fig. 2b) with angles $\alpha$ and $\beta$ and a sketch of the
group velocity together with the asymptote (pink line) of the surface mode
dispersion relation from Eq.(13). The circle with radius $1/(k_{0}z)\approx
21$ determines the location of $(k_{x},k_{y})$ points which give the main
contribution for the inter-particle HF via the surface at a distance of
$z=50$nm.

Figure 5: a) Power $P_{2}$ exchanged between the two NPs (with optical axis of
substrate and NP either in x- or y direction) as function of the rotation
angle $\phi$ of the axis connecting the NPs and the x-axis. The used
parameters are $T_{1}=350$K, $T_{2}=300K$, $z=50$nm, $t=500$nm, $d=700$nm and
$R=10$nm. The power is normalized to the value $P_{2}(\phi_{\rm OA})$ with
$\phi_{\rm OA}=0$ for OA in y-direction and $\phi_{\rm OA}=\pi/2$ for OA in
x-direction. b) Power spectrum $P_{2,\omega}$ for the different rotational
angles with OA in y-direction. Inset: enlargement of the frequency region
around the resonance at $\omega=2.906\times 10^{14}\,{\rm rad/s}$.

Furthermore, the modes which are responsible for the heat transfer when
bringing the NPs in the near-field regime of the slab are those modes for
which $\sqrt{k_{x}^{2}+k_{y}^{2}}\approx 1/z$, because these are the modes
which dominate the near-fields of the slabs Volokitin2007 ; Doro2011 . Hence,
when drawing a circle of radius $1/(k_{0}z)$ into the $q_{x}$-$q_{y}$-plane
then the modes crossing this circles will play the most important role for the
inter-particle HF via the modes of the substrate as shown in Fig. 4 where
$z=50\,{\rm nm}$ and hence $1/(k_{0}z)\approx 21$. It can be seen that there
are significant contributions mainly due to the large wave vector HSPP, but
also by the large wave-vector HVPP which have a crossing with the circle at
the angles $\beta$ or $\alpha$, resp. Hence, when rotating the substrate or
the NPs in the x-y plane one can expect to have a large inter-particle HF when
the inter-particle axis is aligned to the direction of the group velocity

$\mathbf{v}_{g}=\nabla_{\mathbf{k}}\omega=\begin{pmatrix}\frac{{\rm
d}\omega}{{\rm d}k_{x}}\\\ \frac{{\rm d}\omega}{{\rm
d}k_{y}}\end{pmatrix}=\begin{pmatrix}\bigl{(}\frac{{\rm d}k_{x}}{{\rm
d}\omega}\bigr{)}^{-1}\\\ \bigl{(}\frac{{\rm d}k_{y}}{{\rm
d}\omega}\bigr{)}^{-1}\end{pmatrix}$ (12)

of the HSPP which can be evaluated from the dispersion relation of the HSPP.
As depicted in Fig. 4 the group velocity is always perpendicular to the
isofrequency lines. Therefore we can easily determine its direction in the
following approximative manner. In the long-wave vector limit $q_{x}\gg 1$ and
$q_{y}\gg 1$ the HSPP for a film converge to the bulk dispersion relation in
Eq. (9) which can be approximated in this limit by

$\frac{k_{y}}{k_{x}}=\frac{q_{y}}{q_{x}}\approx\sqrt{\frac{\epsilon_{\perp}^{2}-1}{1-\epsilon_{\parallel}\epsilon_{\perp}}}=\tan(\beta)$
(13)

where $\beta$ is the angle between the optical axis and the large wave-vector
branch of the surface mode. If now the axis connecting the NP is along the
y-axis, then we have to rotate the substrate or the NPs by the angle $\beta$
to align the inter-particle axis with the direction of the group velocity. On
the other hand, if the NPs are aligned along the x-axis it is necessary to
rotate the substrate or the NPs by the angle $\alpha$ where
$\tan(\alpha)\tan(\beta)=1$. Since, at $\omega=2.906\times 10^{14}$rad/s we
have ${\rm Re}(\epsilon_{\perp})=-2$ and ${\rm Re}(\epsilon_{\parallel})=2.8$
we find $\beta=0.19\pi$ and $\alpha=0.31\pi$. To verify this reasoning we show
in Fig. 5a) the power $P_{2}$ as a function of the rotation angle $\phi$ of
the substrate when starting with both NPs aligned along the x-axis (and for
two orientations of the optical axis). It can be seen that there is indeed the
maximum at the expected position which is more than 1400 times larger than the
for transmitted power when the axis connecting the NPs is perpendicular to the
OA. Hence, we find a huge modulation contrast due to the strong directionality
of the HSPP. Furthermore, in Fig. 5b) we show that for all rotational angles
the high-frequency resonance dominates the power flow between the NPs in this
configuration. Note, that the direction of the rotation is due to the symmetry
irrelevant.

Figure 6: Normalized spectral power $P_{2,\omega}(\phi)/P_{2,\omega}(\phi=0)$
as a function of the rotation angle $\phi$ for different inter-particle
distances $d$ using parameters as in Fig. (5) and optical axis of substrate
and NPs in y direction. Figure 7: Propagation length $\Lambda$ of the HSPP at
$\omega=2.9\times 10^{14}$rad/s from Eq. (14) for a semi-infinite substrate as
function of $q=\sqrt{q_{x}^{2}+q_{y}^{2}}$. Figure 8: As in Fig. 6 for
$d=700$nm and varying $t$.

In Fig. 6 we show how the modulation of the transferred power $P_{2,\omega}$
at the resonance frequency $\omega=2.9\times 10^{14}$ for different rotational
angles $\phi$ depends on the inter-particle distance $d$ when starting with
NPs axis along x-direction ( $\phi=0$) and optical axis along y direction.
First, it is obvious that there is always a maximum at $\alpha\approx 0.31\pi$
as expected from our explanation. Furthermore, it is clear that for angles
larger than the asymptote of the dispersion relation of the HVPP with the
y-axis also HVPP contribute to the HF so that at $\phi=\pi/2$ the inter-
particle HF is in general larger than for $\phi=0$. It can further be observed
that there is an optimal distance of about $d=1200$nm at which the enhancement
at $\phi=0.31\pi$ has its maximum of about 3500 times $P_{2,\omega}(\phi=0)$.
Note that of course for $d\ll z$ the particle-particle heat transfer is just a
direct heat transfer via the vacuum. By increasing the inter-particle distance
the HVPP and HSPP in the substrate will more and more contribute to the inter-
particle HF until the distance $d$ is getting larger than the propagation
length of the HSPP or HVPP. Then the HF via between the NPs will decay when
further increasing the inter-particle distance.

Figure 9: Imaginary part of $r_{\rm pp}$ for different thicknesses of the
substrate at $\omega=2.9\times 10^{14}$rad/s: a) $t=300$nm, b) $t=200$nm, c)
$t=100$nm, d) $t=50$nm. The dashed lines have the same meaning as in Fig. 2.

To see that the distance $d$ of maximal heat transfer is intimately connected
to the propagation length $\Lambda$ of the HSPP we first define it by

$\Lambda=-\frac{|\mathbf{v}_{\rm g}|}{2\omega^{\prime\prime}}$ (14)

where $\omega^{\prime\prime}$ is the imaginary part of the complex solution of
the dispersion relation in Eq. (9) or (10), resp., assuming real $k_{x}$ and
$k_{y}$. Since the fields of the HSPP have a phase $\exp({\rm i}k_{x}x+{\rm
i}k_{y}y-{\rm i}\omega t)$, $\omega^{\prime\prime}$ is negative and the
damping constant for the intensity is $-1/2\omega^{\prime\prime}$. In Fig. 7
we have plotted the propagation length of the HSPP at $\omega=2.9\times
10^{14}\,{\rm rad/s}$. Obviously, at $q=1/(k_{0}z)\approx 21$ we have a value
of $1.5\,\mu{\rm m}$ which is very good agreement with the maximum found at
$d=1200$nm. Hence, by increasing the distance $z$ of the nanoparticles to the
substrate the maximal HF increase shifts to larger $d$ and vice versa.

Finally, in Fig. 8 we also show the impact of the substrate thickness $t$ on
the spectral transferred power at $\omega=2.9\times 10^{14}$rad/s. It can be
observed that the maximum angle is shifted to smaller angles when making the
substrate thinner. For thicknesses smaller than 300nm clearly a second maximum
at a larger angle can be observed. To understand this behaviour we show in
Fig. 9 the corresponding plots of ${\rm Im}(r_{\rm pp})$ which clearly show
the mode structure available for the inter-particle HF. It can be seen that
the surface and propagating mode dispersion are pushed to larger wave vectors
when making the film thinner. Due to this effect the dispersion of the HSPP in
the first quadrant is getting slightly steeper so that the angle $\beta$ is
slightly increased and hence $\alpha$ is slightly decreased. This explains the
shift of the first maximum in Fig. 9. That the maximum also drops to smaller
value when making the substrate thinner can also be attributed to the fact
that less and less HSPP can contribute since the dispersion relation is pushed
out of the circle with radius $21$. The contributions for angles larger than
the maximum are again due to the HVPP. As for the HSPP for thin layers less
and less of these modes can contribute to the radiative HF explaining the
decreasing inter-particle HF for thicknesses below 300nm and angles larger
than the maximum angle. For very thin substrates of only $t=50$nm these modes
do simply not exist in the wave vector region responsible for the heat
transfer as can be seen in Fig. 9d) explaining the very small HF for angles
larger than the maximum angle for $t=50\,{\rm nm}$. It can also be understood
that there is a sudden HF drop for thin films for angles slightly larger than
the maximum angle. This drop falls into the angle region where we have no HSPP
and no HVPP. When making the films thinner this region becomes larger and can
at a certain thickness directly be observed in the angle dependence of the
inter-particle HF. Hence, we can associate the second maximum to the
contribution of the propagating HVPP. These are of course also the reason for
the maxima at larger angles than $0.31\pi$ in Fig. 6.

##  V conclusion

In summary, we have shown that the radiative HF between two NPs in close
vicinity to a natural hyperbolic material like hBN can be highly modulated by
exploiting their coupling to the highly directional HSPP. The modulation
contrast can be almost 1500 when rotating the axis connecting the NPs with
respect to the optical axis of the substrate. This contrast is especially high
for optically thick films. For thin films with a thickness $t$ on the order of
the distance $z$ or smaller the HF between the NPs cannot couple anymore to
the HSPP so that in this case the HF is efficiently inhibited. Apart from its
potential usefulness for thermal management at the nanoscale the observed
effect might also be employed for a future multitip HF measurement PBA2019 of
the HSPP dispersion relation itself, since it is highly dependent on the
dispersion relation of the HSPP. A direct measurement of the modulation effect
might be possible with a many-body heat transfer setup like that recently
developed by Reddy’s group Reddy .

##  VI Acknowledgment

S.-A. B. acknowledges support from Heisenberg Programme of the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under the project No.
404073166. Xiaohu Wu acknowledges support from Natural Science Foundation of
Shandong Province under the project No. ZR2020LLZ004. Yang Hu acknowledges
support from Key Research and Development Program of Shaanxi Province under
the project No. 2018SF-387 and Basic Research Program of Taicang under the
project No. TC2019JC01.

##  Appendix A Alternative derivation of the angles for maximal heat flux

An alternative route to the foregoing around Eq. (13) is to determine first
the components of the group velocity

$\displaystyle\frac{{\rm d}k_{x}}{{\rm d}\omega}$
$\displaystyle=k_{y}\frac{{\rm d}}{{\rm
d}\omega}\sqrt{\frac{1-\epsilon_{\parallel}\epsilon_{\perp}}{\epsilon_{\perp}^{2}-1}}=k_{y}\frac{{\rm
d}}{{\rm d}\omega}\tan(\alpha),$ (15) $\displaystyle\frac{{\rm d}k_{y}}{{\rm
d}\omega}$ $\displaystyle=k_{x}\frac{{\rm d}}{{\rm
d}\omega}\sqrt{\frac{\epsilon_{\perp}^{2}-1}{1-\epsilon_{\parallel}\epsilon_{\perp}}}=k_{x}\frac{{\rm
d}}{{\rm d}\omega}\tan(\beta)$ (16)

so that the group velocity is

$\mathbf{v}_{g}=\begin{pmatrix}\frac{1}{k_{y}\frac{{\rm d}}{{\rm
d}\omega}\tan(\alpha)}\\\ \frac{1}{k_{x}\frac{{\rm d}}{{\rm
d}\omega}\tan(\beta)}\end{pmatrix}.$ (17)

Therefore, the angle of the group velocity with the y-axis is determined by

$\begin{split}\tan(\tilde{\beta})&=\frac{|k_{x}\frac{{\rm d}}{{\rm
d}\omega}\tan(\beta)|}{|k_{y}\frac{{\rm d}}{{\rm d}\omega}\tan(\alpha)|}\\\
&=\tan(\alpha)\biggl{|}\frac{\frac{{\rm d}}{{\rm
d}\omega}\tan(\beta)}{\frac{{\rm d}}{{\rm d}\omega}\tan(\alpha)}\biggr{|}\\\
&=\tan(\alpha)\biggl{|}-\frac{\frac{{\rm d}\beta}{{\rm
d}\omega}\frac{1}{\sin^{2}(\beta)}}{\frac{{\rm d}\beta}{{\rm
d}\omega}\frac{1}{\cos^{2}(\beta)}}\biggr{|}\\\
&=\tan(\alpha)\tan^{2}(\beta)\\\ &=\tan(\beta).\end{split}$ (18)

Obviously, we obtain the same result as with the graphical method.

##  Appendix B Green Function

The Green function $\mathds{G}_{ij}=\mathds{G}_{ij}^{0}+\mathds{G}_{ij}^{\rm
sc}$ consists of a vacuum part $\mathds{G}_{ij}^{0}$ and a scattering part
$\mathds{G}_{ij}^{\rm sc}$. The vacuum part $\mathds{G}^{0}_{ij}$ is given by
Novotny

$\mathds{G}^{0}_{ij}=\frac{e^{{\rm i}k_{0}{d}}}{4\pi
d}(a\mathds{1}+b\mathbf{e}_{ij}\otimes\mathbf{e}_{ij})$ (19)

with $d=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ and

$a=1+\frac{ik_{0}{d}-1}{k_{0}^{2}d^{2}},b=\frac{3-3ik_{0}{d}-k_{0}^{2}d^{2}}{k_{0}^{2}d^{2}},$
(20)

and

$\mathbf{e}_{ij}=\frac{\mathbf{r}_{i}-\mathbf{r}_{j}}{d}.$ (21)

The scattering part $\mathds{G}_{ij}^{\rm sc}$ is the contribution due to the
presence of the film. Assuming that $z_{i}=z_{j}=z$ it can be expressed within
the Weyl representation as paper_2sic ; Ott2019 ; Ott2020

$\mathds{G}^{sc}_{ij}=\int\frac{{\rm d}^{2}\kappa}{(2\pi)^{2}}e^{{\rm
i}\boldsymbol{\kappa}\cdot(\mathbf{x}_{i}-\mathbf{x}_{j})}\frac{{\rm i}{\rm
e}^{2{\rm
i}\gamma_{0}z}}{2\gamma_{0}}\\!\\!\sum_{\alpha,\beta=p,s}\\!\\!r_{\alpha\beta}\mathbf{a}_{\alpha}^{\,+}\otimes\mathbf{a}_{\beta}^{\,-}$
(22)

with $\gamma_{0}=\sqrt{k_{0}^{2}-\kappa^{2}}$,
$\boldsymbol{\kappa}=(k_{x},k_{y})$, $\mathbf{x}_{i/j}=(x_{i/j},y_{i/j})^{t}$,
and the two polarization vectors

$\mathbf{a}_{s}^{\pm}=\frac{1}{\kappa}\begin{pmatrix}k_{y}\\\ -k_{x}\\\
0\end{pmatrix}\quad\text{and}\quad\mathbf{a}_{p}^{\pm}=\frac{1}{\kappa
k_{0}}\begin{pmatrix}\mp k_{x}\gamma_{0}\\\ \mp k_{y}\gamma_{0}\\\
\kappa^{2}\end{pmatrix}.$ (23)

The double integral over the whole range of $k_{x}$ and $k_{y}$ cannot be
further simplified because due to the anisotropy of the considered material
hBN the reflection coefficients $r_{ss}$, $r_{sp}$, $r_{ps}$, and $r_{pp}$
depend on $\boldsymbol{\kappa}$ and not only $\kappa$ as for isotropic
materials. We evaluate these integrals numerically taking as input the
reflection coefficients of the hBN film which are derived with the method
detailed in App. C. Note, that there is a typo in the expression for the
scattering part of the Green function in Refs. Ott2019 ; Ott2020 : there the
factor 2 in the exponential appears as a prefactor.

##  Appendix C Reflection coefficients

Here, we sketch the derivation of the reflection coefficients for a uni-axial
film with optical axis in a plane parallel to the film interfaces. In the
following the two interfaces of the film of thickness $t$ are assumed to be at
$z=0$ and $z=-t$. Taking the translational symmetry parallel to the film
interfaces into account the electromagnetic field within the film can be
expressed by

$\displaystyle\mathbf{E}_{m}=\mathbf{S}(z)e^{{\rm i}(k_{x}x+k_{y}y-\omega t)}$
(24)
$\displaystyle\mathbf{H}_{m}=\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\mathbf{U}(z)e^{{\rm
i}(k_{x}x+k_{y}y-\omega t)}$ (25)

with $\mathbf{S}=(S_{x},S_{y},S_{z})^{t}$ and
$\mathbf{U}=(U_{x},U_{y},U_{z})^{t}$ which are for the moment unknown. For
further simplification we define $K_{y}\equiv k_{y}/k_{0}$, $K_{x}\equiv
k_{x}/k_{0}$ and $z^{\prime}=zk_{0}$. Now, applying the macroscopic Maxwell
equations

$\displaystyle\nabla\times\mathbf{E}_{m}={\rm i}\omega\mu_{0}\mathbf{H}_{m},$
(26) $\displaystyle\nabla\times\mathbf{H}_{m}=-{\rm
i}\omega\epsilon_{0}\uuline{\epsilon}\mathbf{E}_{m},$ (27)

for a uni-axial medium (with optical axis in x-y plane) described by the
permittivity tensor $\uuline{\epsilon}$, we obtain

$\partial_{z^{\prime}}\begin{pmatrix}S_{x}\\\ S_{y}\\\ U_{x}\\\
U_{y}\end{pmatrix}=\mathds{A}\begin{pmatrix}S_{x}\\\ S_{y}\\\ U_{x}\\\
U_{y}\end{pmatrix}$ (28)

with

${\small\mathds{A}={\rm
i}\begin{pmatrix}0&0&\frac{K_{x}K_{y}}{\epsilon_{zz}}&1-\frac{K_{x}^{2}}{\epsilon_{zz}}\\\
0&0&-1+\frac{K_{y}^{2}}{\epsilon_{zz}}&-\frac{K_{x}K_{y}}{\epsilon_{zz}}\\\
-K_{x}K_{y}-\epsilon_{yx}&K_{x}^{2}-\epsilon_{yy}&0&0\\\
\epsilon_{xx}-K_{y}^{2}&K_{x}K_{y}+\epsilon_{xy}&0&0\end{pmatrix}.}$ (29)

The solution of Eq. (28) can be expressed by linear combinations

$\displaystyle
S_{x}(z)=\sum_{m=1}^{2}w_{1,m}c_{m}e^{k_{0}q_{m}(z-z_{t})}+\sum_{m=1}^{2}w_{1,2+m}c_{m}e^{k_{0}q_{m+2}z}$
(30) $\displaystyle
S_{y}(z)=\sum_{m=1}^{2}w_{2,m}c_{m}e^{k_{0}q_{m}(z-z_{t})}+\sum_{m=1}^{2}w_{2,2+m}c_{m}e^{k_{0}q_{m+2}z}$
(31) $\displaystyle
U_{x}(z)=\sum_{m=1}^{2}w_{3,m}c_{m}e^{k_{0}q_{m}(z-z_{t})}+\sum_{m=1}^{2}w_{3,2+m}c_{m}e^{k_{0}q_{m+2}z}$
(32) $\displaystyle
U_{y}(z)=\sum_{m=1}^{2}w_{4,m}c_{m}e^{k_{0}q_{m}(z-z_{t})}+\sum_{m=1}^{2}w_{4,2+m}c_{m}e^{k_{0}q_{m+2}z}$
(33)

of the four eigenvectors
$\mathbf{w}_{m}=(w_{1,m},w_{2,m},w_{3,m},q_{4,m})^{t}$ ($m=1,2,3,4$) of
$\mathds{A}$ to the four eigenvalues $q_{m}$ with ${\rm Im}(q_{m})>0$
($m=1,2$) and ${\rm Im}(q_{m})<0$ ($m=3,4$). The four coefficients $c_{m}$ are
fixed by imposing the boundary conditions.

To impose the boundary conditions we need also an ansatz for the fields
outside the film. Let’s start with the incident and reflected fields in the
region $z>0$. We make the ansatz

$\displaystyle\mathbf{E}_{i}$
$\displaystyle=\bigl{(}A_{s}\mathbf{a}_{s}^{-}+A_{p}\mathbf{a}_{p}^{-}\bigr{)}{\rm
e}^{{\rm i}(k_{x}x+k_{y}y-k_{z}z-\omega t)},$ (34)
$\displaystyle\mathbf{H}_{i}$
$\displaystyle=\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\bigl{(}-A_{s}\mathbf{a}_{p}^{-}+A_{p}\mathbf{a}_{s}^{-}\bigr{)}{\rm
e}^{{\rm i}(k_{x}x+k_{y}y-k_{z}z-\omega t)}$ (35)

and

$\displaystyle\mathbf{E}_{r}$
$\displaystyle=\bigl{(}R_{s}\mathbf{a}_{s}^{+}+R_{p}\mathbf{a}_{p}^{+}\bigr{)}{\rm
e}^{{\rm i}(k_{x}x+k_{y}y+k_{z}z-\omega t)},$ (36)
$\displaystyle\mathbf{H}_{r}$
$\displaystyle=\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\bigl{(}-R_{s}\mathbf{a}_{p}^{+}+R_{p}\mathbf{a}_{s}^{+}\bigr{)}{\rm
e}^{{\rm i}(k_{x}x+k_{y}y+k_{z}z-\omega t)}$ (37)

with $k_{z}=\sqrt{k_{0}^{2}-\kappa^{2}}$. Similarly, we can make the ansatz
for the transmitted fields in region $z<-t$

$\displaystyle\mathbf{E}_{t}$
$\displaystyle=\bigl{(}T_{s}\mathbf{a}_{s}^{-}+T_{p}\mathbf{a}_{p}^{-}\bigr{)}{\rm
e}^{{\rm i}(k_{x}x+k_{y}y-k_{z}z-\omega t)},$ (38)
$\displaystyle\mathbf{H}_{t}$
$\displaystyle=\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\bigl{(}-T_{s}\mathbf{a}_{p}^{-}+T_{p}\mathbf{a}_{s}^{-}\bigr{)}{\rm
e}^{{\rm i}(k_{x}x+k_{y}y-k_{z}z-\omega t)}.$ (39)

Finally, by demanding that the field components parallel to the interface are
continuous at $z=0$ and $z=-t$ (which are the boundary conditions), we obtain
eight equations which fix the unknown eight coefficients $c_{m}$, $R_{s/p}$,
and $T_{s,p}$ for given $A_{s}$ and $A_{p}$. In particular, the reflection
coefficients are determined by $r_{ss}=R_{s}/A_{s}$, $r_{ps}=R_{p}/A_{s}$ for
an arbitrary $A_{s}\neq 0$ and $A_{p}=0$ and $r_{sp}=R_{s}/A_{p}$,
$r_{pp}=R_{p}/A_{p}$ for an arbitrary $A_{p}\neq 0$ and $A_{s}=0$.

## References

  * (1) K. Sääskilathi, J. Oksanen J. Tulkki, Quantum Langevin equation approach to electromagnetic energy transfer between dielectric bodies in an inhomogeneous environment, Phys. Rev. B 89, 134301 (2014).
  * (2) K. Asheichyk, B. Müller M. Krüger, Heat radiation and transfer for point particles in arbitrary geometries, Phys. Rev. B 96, 155402 (2017).
  * (3) J. Dong, J. Zhao, and L. Liu, Long-distance near-field energy transport via propagating surface waves, Phys. Rev. B 97, 075422 (2018).
  * (4) R. Messina, S.-A. Biehs, and P. Ben-Abdallah, Surface-mode-assisted amplification of radiative heat transfer between nanoparticles, Phys. Rev. B 97, 165437 (2018).
  * (5) M.-J. He, H. Qi, Y.-T. Ren, Y.-J. Zhao, M. Antezza, Giant thermal magnetoresistance driven by graphene magnetoplasmon, Appl. Phys. Lett. 115, 263101 (2019)
  * (6) R. Messina, M. Antezza, P. Ben-Abdallah, Three-body amplification of photon heat tunneling, Phys. Rev. Lett. 109, 244302 (2012).
  * (7) R. Messina, P. Ben-Abdallah, B. Guizal, M. Antezza, S.-A. Biehs, Hyperbolic waveguide for long-distance transport of near-field heat flux, Phys. Rev. B 94, 104301 (2016).
  * (8) Y. Zhang, H.-L. Yi , H.-P. Tan, M. Antezza, Giant resonant radiative heat transfer between nanoparticles, Phys. Rev. B 100, 134305 (2019).
  * (9) S. McSherry and A. Lenert, Extending the thermal near field through compensation in hyperbolic waveguides, Phys. Rev. Applied 14, 014074 (2020).
  * (10) L. Zhu and S. Fan, Persistent directional current at equilibrium in nonreciprocal many-body near field electromagnetic heat transfer, Phys. Rev. Lett. 117, 134303 (2016).
  * (11) L. Zhu, Y. Guo and S. Fan, Theory of many-body radiative heat transfer without the constraint of reciprocity, Phys. Rev. B 97, 094302 (2018).
  * (12) I. Latella and P. Ben-Abdallah, Giant thermal magnetoresistance in plasmonic structures, Phys. Rev. Lett. 118, 173902, (2017).
  * (13) R. M. Abraham Ekeroth, P. Ben-Abdallah, J.C. Cuevas, and A. García-Martín, Anisotropic thermal magnetoresistance for an active control of radiative heat transfer, ACS Photonics 5, 705 (2017).
  * (14) M.-J. He, H. Qi, Y.-X. Su, Y.-T. Ren, Y.-J. Zhao, and M. Antezza, Giant thermal magnetoresistance driven by graphene magnetoplasmon, Appl. Phys. Lett. 117, 113104 (2020).
  * (15) E. Moncada-Villa, V. Fernández-Hurtado, F. J. García-Vidal, A. García-Martín, and J. C. Cuevas, Magnetic field control of near-field radiative heat transfer and the realization of highly tunable hyperbolic thermal emitters, Phys. Rev. B 92, 125418 (2015).
  * (16) J. Song, Q. Cheng, L. Lu, B. Li, K. Zhou, B. Zhang, Z. Luo, and X. Zhou, Magnetically tunable near-field radiative heat transfer in hyperbolic metamaterials, Phys. Rev. Applied 13, 024054 (2020).
  * (17) P. Ben-Abdallah, Photon thermal hall effect, Phys. Rev. Lett. 116, 084301, (2016).
  * (18) A. Ott, S.-A. Biehs, and P. Ben-Abdallah, Anomalous photon thermal Hall effect, Phys. Rev. B 101, 241411(R) (2020).
  * (19) A. Ott, P. Ben-Abdallah, and S.-A. Biehs, Circular heat and momentum flux radiated by magneto-optical nanoparticles, Phys. Rev. B 97, 205414 (2018).
  * (20) C. Khandekar, Z. Jacob, Thermal spin photonics in the near-field of nonreciprocal media, New J. Phys. 21, 103030 (2019).
  * (21) C. Khandekar and Z. Jacob, Circularly polarized thermal radiation from nonequilibrium coupled antennas, Phys. Rev. Appl. 12, 014053 (2019).
  * (22) A. Ott, R. Messina, P. Ben-Abdallah and S.-A. Biehs, Radiative thermal diode driven by nonreciprocal surface waves, Appl. Phys. Lett. 114, 163105 (2019).
  * (23) A. Ott and S.-A. Biehs, Thermal rectification and spin-spin coupling of nonreciprocal localized and surface modes, Phys. Rev. B 101, 155428 (2020).
  * (24) G. Domingues, S. Volz, K. Joulain, and J.-J. Greffet, Heat transfer between two nanoparticles through near field interaction, Phys. Rev. Lett. 94, 085901 (2005).
  * (25) P.-O. Chapuis, M. Laroche, S. Volz, and J.-J. Greffet, Near-field induction heating of metallic nanoparticles due to infrared magnetic dipole contribution, Phys. Rev. B 77, 125402 (2008).
  * (26) Y. Zhang, M. Antezza, H.-L. Yi , H.-P. Tan, Metasurface-mediated anisotropic radiative heat transfer between nanoparticles, Phys. Rev. B 100, 085426 (2019).
  * (27) R. Yu, A. Manjavacas, and F. J. García de Abajo, Ultrafast radiative heat transfer, Nat. Comm. 8, 2 (2017).
  * (28) O. Ilic, N. H. Thomas, T. Christensen, M. C. Sherrott, Marin Soljačić, A. J. Minnich, O. D. Miller, and H. A. Atwater, Active radiative thermal switching with graphene plasmon resonators, ACS Nano 12, 2474 (2018).
  * (29) Y. Huang, S. V. Boriskina, and G. Chen, Electrically tunable near-field radiative heat transfer via ferroelectric materials, Appl. Phys. Lett. 105, 244102 (2014).
  * (30) S.-A. Biehs, R. Messina, P. S. Venkataram, A. W. Rodriguez, J. C. Cuevas, P. Ben-Abdallah, Near-field Radiative Heat Transfer in Many-Body Systems, Rev. Mod. Phys. 93, 025009 (2021).
  * (31) E. E. Narimanov and A. V. Kildishev, Naturally hyperbolic, Nat. Phot. 9, 214 (2015).
  * (32) X. Liu, J. Shen, and Y. Xuan, Near-field thermal radiation of nanopatterned black phosphorene mediated by topological transitions of phosphorene plasmons, Nanosc. Microsc. Therm. Phys. Eng. 23, 188 (2019).
  * (33) J. Shen, X. Liu, and Y. Xuan, Near-field thermal radiation between nanostructures of natural anisotropic material, Phys. Rev. Appl. 10, 034029 (2018).
  * (34) X. Liu and Y. Xuan, Super-Planckian thermal radiation enabled by hyperbolic surface phonon polaritons, Sci China Tech Sci 59, 1680 (2016).
  * (35) X. Wu and C. Fu, Near-field radiative heat transfer between uniaxial hyperbolic media: Role of volume and surface phonon polaritons, JQSRT 258, 107337 (2021).
  * (36) X. Wu and C. Fu, Hyperbolic volume and surface phonon polaritons excited in an ultrathin hyperbolic slab: connection of dispersion and topology, Nanosc. Microsc. Therm. Phys. Eng. 25, 64 (2021).
  * (37) X. Liu, J. Shen, and Y. Xuan, Pattern-free thermal modulator via thermal radiation between Van der Waals materials, JQSRT 200, 100 (2017)
  * (38) X. Wu and Ceji Fu, Int. J., Near-field radiative modulator based on dissimilar hyperbolic materials with in-plane anisotropy, Heat and Mass Transfer 168, 120908 (2021).
  * (39) S.-A. Biehs, F. S. S. Rosa, and P. Ben-Abdallah, Modulation of near-field heat transfer between two gratings, Appl. Phys. Lett. 98, 243102 (2011).
  * (40) R. M. A. Ekeroth, A. Garcia-Martin, and J. C. Cuevas, Thermal discrete dipole approximation for the description of thermal emission and radiative heat transfer of magneto-optical systems, Phys. Rev. B 95, 235428 (2017).
  * (41) A. Lakhtakia, V. K. Varadan, and V. V. Varadan, Low-frequency scattering by an imperfectly conducting sphere immersed in a dc magnetic field, International Journal of Infrarared and Millimeter Waves 12, pp. 1253-1264 (1991).
  * (42) P. Li, I. Dolado, F. J. A. Mozaz, A. Y. Nikitin, F. Casanova, L. E. Hueso, S. Velez, and R. Hillenbrand, Optical nanoimaging of hyperbolic surface polaritons at the edges of van der Waals materials, Nano Lett. 17, 228-235 (2017).
  * (43) X. H. Wu, C. J. Fu, and Z. M. Zhang, Influence of hBN orientation on the near-field radiative heat transfer between graphene/hBN heterostructures, J. Phot. En. 9, 032702 (2018).
  * (44) G. Alvarez-Perez, K. V. Voronin, V. S. Vokov, Analytical approximations for the dispersion of electromagnetic modes in slabs of biaxial crystals, Phys. Rev. B 100, 235408 (2019).
  * (45) A. I. Volokitin and B. N. J. Persson, Near-field radiative heat transfer and noncontact friction, Rev. Mod. Phys. 79, 1291 (2007).
  * (46) I. A. Dorofeyev and E. A. Vinogradov,Fluctuating electromagnetic fields of solids, Phys. Rep. 504 75 (2011).
  * (47) P. Ben-Abdallah, Multitip near-field scanning thermal microscopy, Phys. Rev. Lett. 123, 264301 (2019).
  * (48) D. Thompson, L. Zhu, E. Meyhofer, and P. Reddy, Nanoscale radiative thermal switching via multi-body effects, Nat. Nanotechn. 15, 99 (2020).
  * (49) L. Novotny and B. Hecht, Principles of Nano-Optics, Cambridge University Press, (2006).


* Huix-Rotllant et al. 2010 Huix-Rotllant, M.; Natarajan, B.; Ipatov, A.; Wawire, C. M.; Deutsch, T.; Casida, M. E. Assessment of noncollinear spin-flip Tamm–Dancoff approximation time-dependent density-functional theory for the photochemical ring-opening of oxirane. _Physical Chemistry Chemical Physics_ 2010, _12_ , 12811–12825
  * Neese et al. 2009 Neese, F.; Hansen, A.; Liakos, D. G. Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis. _The Journal of chemical physics_ 2009, _131_ , 064103
  * Neese et al. 2009 Neese, F.; Hansen, A.; Wennmohs, F.; Grimme, S. Accurate theoretical chemistry with coupled pair models. _Accounts of chemical research_ 2009, _42_ , 641–648
  * Neese et al. 2009 Neese, F.; Wennmohs, F.; Hansen, A. Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method. _The Journal of chemical physics_ 2009, _130_ , 114108
  * Liakos et al. 2011 Liakos, D. G.; Hansen, A.; Neese, F. Weak molecular interactions studied with parallel implementations of the local pair natural orbital coupled pair and coupled cluster methods. _Journal of Chemical Theory and Computation_ 2011, _7_ , 76–87
  * Hansen et al. 2011 Hansen, A.; Liakos, D. G.; Neese, F. Efficient and accurate local single reference correlation methods for high-spin open-shell molecules using pair natural orbitals. _The Journal of chemical physics_ 2011, _135_ , 214102
  * Riplinger and Neese 2013 Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. _The Journal of chemical physics_ 2013, _138_ , 034106
  * Riplinger et al. 2013 Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural triple excitations in local coupled cluster calculations with pair natural orbitals. _The Journal of chemical physics_ 2013, _139_ , 134101
  * Riplinger et al. 2016 Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse maps–A systematic infrastructure for reduced-scaling electronic structure methods. II. Linear scaling domain based pair natural orbital coupled cluster theory. _The Journal of chemical physics_ 2016, _144_ , 024109
  * Datta et al. 2016 Datta, D.; Kossmann, S.; Neese, F. Analytic energy derivatives for the calculation of the first-order molecular properties using the domain-based local pair-natural orbital coupled-cluster theory. _The Journal of Chemical Physics_ 2016, _145_ , 114101
  * Saitow et al. 2017 Saitow, M.; Becker, U.; Riplinger, C.; Valeev, E. F.; Neese, F. A new near-linear scaling, efficient and accurate, open-shell domain-based local pair natural orbital coupled cluster singles and doubles theory. _The Journal of chemical physics_ 2017, _146_ , 164105
  * Becke 1993 Becke, A. D. Density‐functional thermochemistry. III. The role of exact exchange. _The Journal of chemical physics_ 1993, _98_ , 5648
  * Dunning Jr 1989 Dunning Jr, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. _The Journal of chemical physics_ 1989, _90_ , 1007–1023
  * Pavošević et al. 2017 Pavošević, F.; Peng, C.; Pinski, P.; Riplinger, C.; Neese, F.; Valeev, E. F. SparseMaps–A systematic infrastructure for reduced scaling electronic structure methods. V. Linear scaling explicitly correlated coupled-cluster method with pair natural orbitals. _The Journal of chemical physics_ 2017, _146_ , 174108
  * Hehre et al. 1986 Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. Ab initio Molecular Orbital Theory. 1986
  * Whitten 1973 Whitten, J. L. Coulombic potential energy integrals and approximations. _The Journal of Chemical Physics_ 1973, _58_ , 4496–4501
  * Baerends et al. 1973 Baerends, E.; Ellis, D.; Ros, P. Self-consistent molecular Hartree-Fock-Slater calculations I. The computational procedure. _Chemical Physics_ 1973, _2_ , 41–51
  * Dunlap et al. 1979 Dunlap, B. I.; Connolly, J.; Sabin, J. On some approximations in applications of X$\alpha$ theory. _The Journal of Chemical Physics_ 1979, _71_ , 3396–3402
  * Van Alsenoy 1988 Van Alsenoy, C. Ab initio calculations on large molecules: The multiplicative integral approximation. _Journal of computational chemistry_ 1988, _9_ , 620–626
  * Kendall and Früchtl 1997 Kendall, R. A.; Früchtl, H. A. The impact of the resolution of the identity approximate integral method on modern ab initio algorithm development. _Theoretical Chemistry Accounts_ 1997, _97_ , 158–163
  * Eichkorn et al. 1995 Eichkorn, K.; Treutler, O.; Öhm, H.; Häser, M.; Ahlrichs, R. Auxiliary basis sets to approximate Coulomb potentials. _Chemical physics letters_ 1995, _240_ , 283–290
  * Eichkorn et al. 1997 Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Auxiliary basis sets for main row atoms and transition metals and their use to approximate Coulomb potentials. _Theoretical Chemistry Accounts_ 1997, _97_ , 119–124
  * Neese et al. 2009 Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Efficient, approximate and parallel Hartree-Fock and hybrid DFT calculations. A ‘chain-of-spheres’ algorithm for the Hartree–Fock exchange. _Chemical Physics_ 2009, _356_ , 98–109
  * Lee and Taylor 1989 Lee, T. J.; Taylor, P. R. A diagnostic for determining the quality of single-reference electron correlation methods. _International Journal of Quantum Chemistry_ 1989, _36_ , 199–207
  * Gozem et al. 2013 Gozem, S.; Krylov, A. I.; Olivucci, M. Conical intersection and potential energy surface features of a model retinal chromophore: Comparison of EOM-CC and multireference methods. _Journal of chemical theory and computation_ 2013, _9_ , 284–292