# HUE: Pretrained Model and Dataset for
Understanding Hanja Documents of Ancient Korea

Haneul Yoo1, Jiho Jin1, Juhee Son1, JinYeong Bak2, Kyunghyun Cho3, Alice Oh1
1KAIST, South Korea, 2Sungkyunkwan University, South Korea, 3New York
University, USA
{haneul.yoo, jinjh0123<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Historical records in Korea before the 20 century were primarily written in
Hanja, an extinct language based on Chinese characters and not understood by
modern Korean or Chinese speakers. Historians with expertise in this time
period have been analyzing the documents, but that process is very difficult
and time-consuming, and language models would significantly speed up the
process. Toward building and evaluating language models for Hanja, we release
the Hanja Understanding Evaluation dataset consisting of chronological
attribution, topic classification, named entity recognition, and summary
retrieval tasks. We also present BERT-based models continued pretraining on
the two major corpora from the 14 to the 19 centuries: the Annals of the
Joseon Dynasty and Diaries of the Royal Secretariats. 111All codes, models,
and dataset are available at https://github.com/haneul-yoo/HUE.git We compare
the models with several baselines on all tasks and show there are significant
improvements gained by training on the two corpora. Additionally, we run zero-
shot experiments on the Daily Records of the Royal Court and Important
Officials (DRRI). The DRRI dataset has not been studied much by the
historians, and not at all by the NLP community.

##  1 Introduction

Language | Sentence
---|---
Hanja |  上御經筵。
Modern Korean |  임금이 경연에
|  나아갔다.
Simplified Chinese |  国王参加了皇家讲座。
Traditional Chinese |  國王參加了皇家講座。
English | The King attended
| the Royal Lecture.
Table 1: Example sentence in AJD

Large-scale historical records in Korea were mostly produced during the Joseon
dynasty (1392-1897), and the Institute for the Translation of Korean Classics
(ITKC) keeps a comprehensive database of Korean classics at a scale of
approximately 9 billion characters. This digital archive is a great resource
for Korean historians, but the documents remain in the original Hanja language
222Hanja is a set of characters (script) used in ancient Korean, while Hanmun
is a writing style (language) in the same era. However, we will refer to Hanja
as a language following the conventions of the previous works.. Hanja is an
extinct language, and as Table 1 illustrates with a simple sentence, Hanja is
lexically and syntactically different from modern Korean, as well as
simplified and traditional Chinese. Understanding the documents in the digital
archive is thus difficult and would benefit greatly from a Hanja language
model which can also be used to accelerate the expert translation Vale de Gato
(2015). There are two corpora we can use to train the language model, the
Annals of Joseon Dynasty (AJD), first introduced to the NLP community in (Bak
and Oh, 2015), and the Diaries of the Royal Secretariats (DRS) Kang et al.
(2021).

In this paper, we provide the HUE (Hanja Understanding Evaluation) dataset
consisting of chronological attribution, topic classification, named entity
recognition and summary retrieval, a suite of tasks to help build and evaluate
the Hanja language model. In addition to AJD and DRS, we also work with the
Daily Records of the Royal Court and Important Officials (DRRI). Unlike AJD
and DRS that have been analyzed by historians and contain their annotations,
DRRI lacks such systematic analysis, and we use it for zero-shot learning and
introduce it to the NLP community.

We also provide pretrained language models (PLMs) for Hanja trained on AJD and
DRS, fine-tuned for each task in HUE. Our pretrained models on the corpora
from that era outperform the existing language models built for ancient
Chinese, confirming the need for specially-trained Hanja language models. We
also run additional experiments based on the analyses of entity- and word-
level changes on AJD by controlling input conditions by masking named entity
and giving the time period as input. Finally, we demonstrate the effectiveness
of our Hanja language model for analyzing unseen documents, running zero-shot
experiments for chronological attribution and named entity recognition tasks
on DRRI.

Our main contributions are as follows:

  * •

We release the HUE dataset and Hanja PLMs to support historians to understand
and analyze a large volume of historical documents written in Hanja. To the
best of our knowledge, this is the first work proposing Hanja language models
and releasing a NLP benchmark dataset for ancient Hanja documents.

  * •

We demonstrate that providing key information such as named entity and
document age as input improves the performance of Hanja language model on the
HUE tasks.

  * •

We run zero-shot experiments on several HUE tasks from DRRI which have not
been discussed in the NLP community, and demonstrate the performance of our
Hanja language models on unseen historical documents.

##  2 Background

###  2.1 Hanja

Hanja, the writing system based on ancient Chinese characters, was the main
writing system in Korea before the 20 century, while Hangul, the unique Korean
alphabet, has been the main writing system in Korea from the last century.
Formal records from the Joseon dynasty (1392-1897) are written in Hanja, while
spoken language and some written documents were in Hangul, developed in the
15th century. This co-existence of the written and colloquial languages has
led Hanja to evolve to have the basic syntax of classical Chinese, mixing with
the lexical, semantic, and syntactic characteristics of colloquial Korean.

Hanja is significantly different from both modern Korean and modern Chinese.
Modern Korean uses a different alphabet and structure, and traditional Chinese
shares some characters with Hanja, while the lexicon has evolved greatly to
reflect the temporal, geographical, and cultural differences between the
Joseon dynasty and modern-day China. Simplified Chinese, the current written
language in China has diverged more because of the simplification of the
characters. These differences between Hanja and other related languages would
lead to suboptimal performance when adopting the current Chinese language
models to NLU tasks for the Korean historical Hanja documents.

###  2.2 Dataset

Dataset | Size | Training data | Downstream Tasks | Zero-shot | King
---|---|---|---|---|---
AJD | 230K | ✔ | CA, TC, NER | - | Taejo (1) - Sunjong (27)
DRS | 1,380K | ✔ | - | - | Injo (16) - Sunjong (27)
DRRI | 426K | - | SR | CA, NER | Yeongjo (21) - Sunjong (27)
Table 2: Source corpora chosen for building HUE dataset and PLMs

We describe the three corpora of records written in Hanja during the Joseon
dynasty, whose contents and additional information such as topic and named
entities are provided by historians in IKTC 333https://db.itkc.or.kr/. Table 2
shows the list of the Hanja corpora used.

##### Annals of the Joseon Dynasty (AJD)

also called Veritable Records of the Joseon Dynasty, is a corpus of 27 sets of
chronological records, and each set covers one ruler’s reign. AJD has been
translated into Korean from 1968 to 1993 and includes relevant tags such as
the named entities and dates of the documents
444http://esillok.history.go.kr/. We use AJD for both training our Hanja
language models and building the HUE dataset of NLP tasks.

##### Diaries of the Royal Secretariat (DRS)

is a corpus of detailed records of daily events and official schedules of the
court from the first King Taejo to the last (27) Sunjong. Many of the earlier
records were lost, and we use the extant records starting from the 16 King
Injo. DRS is known to hold the largest amount of authentic historic records
and state secrets of the Joseon Dynasty 555http://english.cha.go.kr/. We use
DRS to continue pretraining the language models.

##### Daily Records of the Royal Court and Important Officials (DRRI)

is a corpus of journals written from the 21 King Yeongjo to the last Emperor
Sungjong and presumably initiated from the diaries of the crown prince who
became the 22 King Jeongjo after he was enthroned. DRRI has official daily
records from both the central and the local governments, so encompasses all
events in the country and reports to the king with summaries. DRRI is known to
include details and events of the late Joseon Dynasty not recorded in the AJD
or DRS 5, thus making it a good corpus for zero-shot experiments. We use DRRI
for the supervised summary retrieval task and the zero-shot experiments for
chronological attribution and named entity recognition.

##  3 HUE Dataset

The HUE (Hanja Understanding Evaluation) dataset is built to assist history
scholars to understand Korean historical records written in Hanja. HUE
consists of chronological attribution, topic classification, named entity
recognition, and summary retrieval, which are tasks that can provide helpful
information for studying the documents. We expect that the language models
based on HUE will ultimately help historians to interpret unseen historical
documents and public to grasp basic concept of those documents. We describe
each task in detail below.

###  3.1 Task Description

##### Chronological Attribution (CA)

is a classification task predicting the ruling king when the document was
written. When given a Hanja document from AJD, a classifier outputs one of the
27 kings of the Joseon dynasty. Chronological attribution of the undiscovered
document is the first step in anthology to interpret and translate it. Korean
historians mostly divide the history of the Joseon Dynasty based on the
reigning king, so that we treat chronological attribution as a classification
task.

##### Topic Classification (TC)

is a multi-class and multi-label classification task to find the topics of the
given document. For TC, we use Hanja document from AJD. We suggest two levels
of topics, namely major and minor categories. The major categories consist of
4 broad topics: politics, economy, society, and culture. The minor categories
go with 106 sub-topics such as diplomacy, agriculture, and science.

##### Named Entity Recognition (NER)

is a sequence tagging task, identifying the two types of named entities,
person and location, from the Hanja document from AJD. We divide train,
validation, and test sets such that there are no overlapping entities across
the sets.

##### Summary Retrieval (SR)

is a task to find the most relevant summary that matches the content among the
summary candidates. For this task, we use DRRI, in which each document is a
pair of summary (_gang_) and detailed content (_mok_). Among 426k articles,
265k articles in DRRI dataset contain both _gamg_ and _mok_. Also, we exclude
those with _gang_ longer than _mok_ , in which _gang_ is not the summary of
_mok_. The final dataset contains 213K pairs of content and the corresponding
summaries. We describe more details of the preprocessing in the Appendix.

##  4 Hanja Pretrained Model

As far as we know, there have been no pretrained language models for the Hanja
language. One can use related LMs, the pretrained models for ancient Chinese
as well as multilingual BERT Devlin et al. (2019) which includes traditional
Chinese in its training corpus. AnchiBERT Tian et al. (2021) is pretrained in
ancient Chinese with the Chinese anthologies written around 1000BC to 200BC.
There is some vocabulary overlap between the Hanja documents and traditional
Chinese corpora, we can adopt multilingual BERT and AnchiBERT to learn the
representations of the Hanja texts.

We propose the pretrained language models suitable for Hanja documents by
continuing pretraining those two models on both AJD and DRS. Table 3 shows the
ratio of unknown tokens in the test set of AJD by each model. It implies that
existing AnchiBERT and multilingual BERT can also be exploited as language
models for Hanja documents written in the Joseon dynasty, but the second phase
of pretraining on the corpora of that era remarkably decreases the ratio of
unknown tokens.

| AnchiBERT | mBERT
---|---|---
| Tian et al. (2021) | Devlin et al. (2019)
original | 0.88% | 0.76%
+AJD/DRS | 0.04% | 0.04%

Table 3:  Unknown token ratio of each model in test set for CA, TC, and NER
task in HUE. The first row indicates the results of the original PLMs without
any additional pretraining, and the second row indicates those with continued
pretraining on AJD and DRS.

##  5 Experiment

###  5.1 Experimental Settings

We conduct experiments on HUE with our pretrained model. For the baseline
model, we use BERT without pretraining and compare it to various BERT models
described in Section 3.1. Specifically, we fine-tune each model to act as a
re-ranker in the retrieval task for the summary retrieval. We first retrieve
top-$k$ relevant _gang_ s (summaries) with BM25 Robertson and Zaragoza (2009)
among all _gang_ s in the dataset. Then, we fine-tune the model as a binary
classifier to determine whether the summary matches the content of the _mok_
with the cross-entropy loss Nogueira and Cho (2019). If the ground truth
summary is not included in the top-$k$ relevant summary candidates, we replace
the last $k^{\text{th}}$ summary with the ground truth. We use $k=12$ for
training and $k=100$ for inference. As the representative metrics, we present
F1 scores for CA, TC, and NER tasks, and Mean Reciprocal Rank (MRR) for SR.
The detailed results with other metrics are also available in Appendix.

###  5.2 Overall Results

| CA | TC | NER | SR
---|---|---|---|---
|  | Major | Minor | Person | Location |
| F1 | F1 | F1 | F1 | F1 | mrr
BERT not pretrained | 54.26 | 68.91 | 61.52 | 92.13 | 87.10 | 52.85
mBERT Devlin et al. (2019) | 75.29 | 79.59 | 76.46 | 91.63 | 86.02 | 67.06
AnchiBERT Tian et al. (2021) | 75.74 | 85.81 | 75.22 | 93.28 | 88.01 | 67.92
mBERT + AJD/DRS | 77.77 | 87.13 | 77.84 | 92.83 | 87.90 | 73.88
AnchiBERT + AJD/DRS | 79.33 | 88.33 | 78.10 | 93.13 | 87.91 | 74.29
Table 4: Evaluation results of PLMs on HUE dataset

Table 4 shows the overall experimental results. The bold and the underlined
texts in the table specify the best and the second best result, respectively.
BERT without any pretraining shows the poorest results across all the tasks.
AnchiBERT and mBERT, which are existing language models on the relevant
domains, show better results, and the models continued pretraining on Hanja
documents achieve the best performance among all tasks. This tendency
indicates that all these tasks on understanding historical documents require
pretraining language models on time-specific and domain-specific data.

AnchiBERT pretrained on the Hanja corpora shows slightly better performance
than mBERT pretrained on the same corpora. We assume this is because the
original training corpora of AnchiBERT are much closer to the Hanja documents,
even the era of those two corpora are completely different. The writing style
of both Hanja documents in the Joseon dynasty and anthologies in ancient China
had come from Classical Chinese and share similarities. On the other hand, the
training corpora of mBERT is a contemporary texts whose characters contains
Traditional Chinese, but the structure and the format might be considerably
changed.

##### Chronological Attribution (CA)

Our models continued pretraining on Hanja corpora outperform other baselines
on CA. Detailed analysis on CA result is illustrated in Section 6.1.

##### Topic Classification (TC)

() Div 1 (4 classes) (a) Div 2 (106 classes)

Figure 1: ROC Curve and AUC for Topic Classification. Each value in the legend
indicates the AUC score.

Figure 1 gives ROC curves and AUC values of each model on each task. Our
models show the similar trends to the overall results, outperforming other
language models. For the evaluation results including F1 score, we find and
set the best threshold to each label by Youden’s index (Youden, 1950).

While F1 score goes down as the number of classes increases from 4 to 106,
there is no significant difference on AUC value. This might result from
consistently high recall achieving around 90% on both tasks. It indicates that
the threshold is too low and models tend to predict plausible topics as many
as they can, which might be solved by controlling the threshold. AnchiBERT
pretrained on AJD and DRS, which shows the best performance, predicts 6.39
labels in average, while the average number of ground truth labels of the
minor categories is 1.97. This is probably due to the meaning overlaps in
minor categories. For instance, minor categories such as revenue, finance,
general price level, and commerce are the sub-categories of economy in the
major categories, whose use case cannot be strictly distinguished. It would be
more appropriate in this case to provide all plausible topics roughly rather
than suggesting the one only with high certainty.

##### Named Entity Recognition (NER)

NER also indicates similar trends to the overall benchmark tasks, but with a
small gap among models including BERT without pretraining. It implies that NER
in Hanja documents is a comparably easy task. This might result from certain
patterns in named entities in Hanja. Most of the person entities are 3 letters
starting with the common characters (family name), and most of the location
entities end with the common characters meaning locations or buildings. All
models tend to predict person entities better than location entities.

##### Summary Retrieval (SR)

All fine-tuned BERT-based re-rankers outperform BM25 whose MRR is merely
29.87%, mostly retrieving the ground truth answer at the first trial.
Likewise, our models shows the best results, while BERT without pretraining
show the lowest MRR. It additionally implies that BERT-based re-rankers might
be exploited for retrieving relevant documents from different chronicles in
terms of written style or contents.

###  5.3 Effect of Entity and Document Age on Language Model

We investigate whether providing additional information as input can improve
the performances of language models. For CA, we mask all named entities in the
input and fine-tuned the language model on masked data to examine the impact
of entity information. For TC and NER, we concatenate document age information
to the input text and run comparative experiments to verify the importance of
time period on historical texts.

##### Entity-Masked Chronological Attribution

| CA (Masked)
---|---
BERT not pretrained | 75.07 ($\Delta$ 20.81)
mBERT Devlin et al. (2019) | 83.44 ($\Delta$ 8.15)
AnchiBERT Tian et al. (2021) | 82.45 ($\Delta$ 6.71)
mBERT + AJD/DRS |  83.57 ($\Delta$ 5.80)
AnchiBERT + AJD/DRS |  83.58 ($\Delta$ 4.25)
Table 5: F1 scores on Chronological Attribution given named entities masked.
The value inside parenthesis indicates the increase in performance after
masking the named entities.

Table 5 shows the difference on experimental results in CA when the named
entities in the given input texts are masked. Compared to the default settings
which does not mask named entities, all models show significant improvements.
This is probably because models can truly focus on the content and changes in
writing style without any disturbance of location entities consistently used
for the whole era. Fine-tuned models with entity masked inputs also achieve
nearly the same level of performance in the inference with plain inputs
including named entities. It suggest to fine-tune models masking named entity
in CA, considering the real scenario whose inference texts lack time period
information.

##### Topic Classifciation and Named Entity Recognition with the Age of
Document

| TC | NER
---|---|---
| Major | Minor | Person | Location
BERT not pretrained | 86.99 ($\Delta$ 18.08) | 70.58 ($\Delta$ 9.29) | 91.90 ($\Delta$ -0.23) | 86.43 ($\Delta$ -0.68)
mBERT Devlin et al. (2019) | 93.57 ($\Delta$ 13.98) | 73.64 ($\Delta$ 3.63) | 93.71 ($\Delta$ 2.08) | 87.84 ($\Delta$ 1.83)
AnchiBERT Tian et al. (2021) | 89.32 ($\Delta$ 3.51) | 73.57 ($\Delta$ 4.27) |  94.82 ($\Delta$ 1.54) |  89.84 ($\Delta$ 1.83)
mBERT + AJD/DRS | 90.02 ($\Delta$ 2.89) | 74.84 ($\Delta$ 3.56) | 94.88 ($\Delta$ 2.05) | 89.88 ($\Delta$ 1.98)
AnchiBERT + AJD/DRS |  90.02 ($\Delta$ 1.69) |  74.54 ($\Delta$ 2.47) | 94.70 ($\Delta$ 1.57) | 89.46 ($\Delta$ 1.54)
Table 6: F1 scores on topic classification and named entity recognition given
document age. The value inside parenthesis indicates the difference on
performance after providing document age.

It is for granted to regard that historical texts written over several eras
reveal the time changes with respect to lexical choices and contents. Section
6.2 confirms the hypothesis above in terms of $n$-grams changes over time.
Table 6 shows gaps between experimental results of TC and NER given which king
reigned when the document was written. Providing document age definitely
increases the performance of classifying topics and tagging named entities.
There was a big gap on difference with non-pretrained BERT in TC, which is
probably due to the poor performance of itself in the original setting. All
models show similar trends on both tasks with improved performances compared
to the original settings without document age as input. It is an obvious
result considering that the first step for ancient manuscript is assuming the
written era. It implies the significance of chronological attribution task in
HUE, conveying that chronological attribution task might improve the
performance of other HUE tasks.

###  5.4 Zero-shot Experiment

| CA | NER
---|---|---
|  | Person | Location
BERT not pretrained | 26.94 | 68.50 | 39.86
mBERT Devlin et al. (2019) | 32.19 | 72.08 | 83.36
AnchiBERT Tian et al. (2021) | 30.65 | 71.85 | 77.23
mBERT + AJD/DRS | 35.28 | 88.48 | 72.08
AnchiBERT + AJD/DRS | 35.85 | 82.86 | 76.53
Table 7: F1 scores of chronological attribution and named entity recognition
task on DRRI in zero-shot settings

Countless number of Hanja documents still remain without any analysis and new
documents continue to be unearthed. Therefore, we run zero-shot experiments to
verify the effectiveness of our language models on extracting information from
the historical documents irrelevant to the training corpora. We use DRRI
dataset which is not included in both pretraining and fine-tuning data of our
Hanja language models and execute CA and NER.

Table 7 shows experimental results with CA and NER on DRRI. All models perform
comparably well on the both tasks, regarding that random model will achieve
approximately 3.70% performances with 27 classes in CA. Also, all models in CA
commonly show high precision which might be due to the monotonous and
redundant phrases in the veritable records. It shows similar trends compared
to the Table 4, but the gap among models was notably emphasized in the zero-
shot settings. This results imply that our CA models might be exploited for
the time period prediction of unseen documents in anthology with a reliable
level.

Our models outperform others on NER achieving absolutely high performances,
though entity maps between AJD and DRRI do not match strictly. Interestingly,
all models tend to predict location entities better than person entities which
is the opposite result compared to the original NER on AJD. It is probably due
to the characteristics of each entity, where location entities are all
commonly used in nationwide while person entity might differ by situation.
Further analysis on person and location entities in the view of time changes
is described in Section 6.2. We present that our models trained on the corpora
of the Joseon dynasty provide reliable results on unseen records, implying
that our model can be exploited for the low-resourced documents.

##  6 Further Analysis

###  6.1 Do Historical Events Affect Language Models?

(a) Log-scale Confusion Matrix (b) Mean Absolute Error per King

Figure 2: Chronological attribution results with AnchiBERT + AJD/DRS. Figure
(a) shows the confusion matrix. In Figure (b), each bar indicates the mean
absolute error and the line indicates the number of samples for each king.

To figure out the effect of historical events on models’ prediction, we
analyzed the output of language models on CA. Figure 2 (a) shows a log-scale
confusion matrix of AnchiBERT continued pretraining on AJD and DRS, and Figure
2 (b) indicates the mean absolute error between the predicted king order and
the ground truth per each King. The x-axis in the Figure 2 (b) means the
changes of time by the king reign period. Each bar in Figure 2 (b) indicates
the mean absolute error between the order of ground truth king and the one of
predicted, and the line graph means the number of samples in the test set.

The results of the last two Emperors, Gojong and Sunjong, are remarkable in
that the model rarely gets confused with those two labels to others and tends
not to fail, showing notable difference on Figure 2 (a) and significantly low
mean absolute error on Figure 2 (b). We believe that this is because our model
learns the difference between those two records and the others in the
historical view and get cues to distinguish them. The last two records are not
treated as AJD in general, since those records were inspected and produced by
the Office of Governer-General during the Japanese colonial period with the
view of Empire of Japan who ruled Korea 666http://esillok.history.go.kr/.

The mean absolute error of the predicted order of each king achieved the
difference around one, except for the first King Taejo and the second King
Jeongjong, whose errors are almost doubled. We hypothesize that this is mainly
because there are too small number of examples in those classes. A similar
tendency where the more samples are, the less mean absolute error be has been
observed in other classes. Also, the writing style of AJD had settled down
from the third King Taejong, and those noisy records might confuse models to
predict the exact dates.

###  6.2 Do Time Changes Affect Written Texts?

It stands to reason that AJD written in five centuries reveal the features of
the language changes. In this section, we investigate the hypothesis above
with respect to the named entities and $n$-grams.

##### Words Change over Time

(a) Taejo (1 King) (b) Seongjong (9 King) (c) Hyojong (17 King) (d) Cheoljong
(25 King)

Figure 3: Number of overlapped trigrams per king era. It shows the changes of
trigrams over kings whose x-axis shows the kings and the y-axis shows the
number of overlapped trigrams.

We analyze how frequently words change over time. For each king, we plot how
many trigrams overlap by each king era in the order. Figure 3 shows overlapped
trigrams in the 1st, 9th, 17th, and 25th king and the detailed results with
all kings are described in Appendix. It is consistently observed that the
closer the king era is, the more trigrams are overlapped. These changes result
from not named entity but lexical choices, considering that person and
location entities account for 6.38% and 2.05% in the characters of AJD,
respectively. It verifies that words used in Joseon dynasty had changed over
time gradually, and it enables the language models to capture those features.

##### Named Entity Changes over Time

(a) Person (b) Location

Figure 4: Relative frequency change of top-10 named entities per king. Each
line indicates the change of the relative frequency of one entity over time,
and the color of the line indicates the king era in which the entity is
contained in the top-10 entities. The x-axis represents the time (the kings in
the order) and the y-axis represents the relative frequency in each king era.

We investigate how the named entities had been used over time. In particular,
we show frequency rates of top-10 frequently-used named entities by each king
era and how they change over time in Figure 4. It implies a strong correlation
between person entity and the passage of time, while there is no explicit
correlation to location entity. Most person entities include officials of the
time or the previous kings, relevant to the time. In contrast, most location
entities include neighboring countries or place names in the Joseon, which are
less dependent on the time. The examples of frequently appeared named entities
are described in Appendix.

##  7 Related Work

ML based NLP techniques have been recently applied to anthology to discover
historical documents such as authorship attribution (Ouamour and Sayoud, 2012;
Sayoud and Ouamour, 2017; Reisi and Mahboob Farimani, 2020; Hossain et al.,
2020), NER (Won et al., 2018; Palladino et al., 2020), and manuscript age
detection (Adam et al., 2018). Assael et al. (2022) proposed Ithaca to restore
ancient Greek inscriptions and perform geographical attribution and
chronological attribution of them.

Along with these works, several works provide language models suited for
historical texts in ancient languages and evaluate those models on existing
NLU tasks, which aims to support understanding those documents considering
that the target languages are mostly extinct. Bamman and Burns (2020) propose
Latin BERT for part-of-speech tagging in ancient Latin script. Tian et al.
(2021) suggest AnchiBERT and evaluate their model on some NLP tasks including
poem topic classification.

However, there has been no research attempting to propose language models in
Hanja, which is a dead language in Korea but absolutely necessary to explore
Korean history. Most of the studies with Hanja only shed lights on translating
historical Hanja documents and use AJD as their corpus (Park et al., 2020; Jin
et al., 2020; Kang et al., 2021).

##  8 Conclusion

We present HUE (Hanja Understanding Evaluation) dataset and BERT-based
pretrained language models for classical Hanja documents. HUE dataset includes
diverse tasks that can support analyzing historical documents written in Hanja
which is an extinct language in Korea: Chronological Attribution (CA), Topic
Classification (TC), Named Entity Recognition (NER), and Summary Retrieval
(SR). Our models pretrained on Hanja corpora outperform other language models
and we observe their performance on zero-shot settings with DRRI which is the
dataset never been introduced in NLP community. The experimental results in
king prediction imply that our models capture the historical events or facts
disclosed in the texts. We also explore several methods to support Hanja
language models such as masking named entities and giving document age as
input sources, based on the analyses on textual features in AJD.

Help of adequate resources in Hanja documents might could fill some caveats in
our work which lacks additional experiments and analyses on the records of
different genre such as poetry, novel, and humanities resulting from the low
resources that we can exploit. However, we expect that our dataset and
accompanying language models might facilitate future works on historical
documents written in Hanja by providing fundamental resources to leverage
unknown Hanja corpora.

## Acknowledgements

We would like to thank Yoonman Heo (Institute for the Translation of Korean
Classics) providing expertise on hanja and Korean Classical Chinese. This
research was supported by the Engineering Research Center Program through the
National Research Foundation of Korea (NRF) funded by the Korean Government
MSIT (NRF-2018R1A5A1059921). This work was partly supported by Institute of
Information & communications Technology Planning & Evaluation (IITP) grant
funded by the Korea government (MSIT) (No.2019-0-00421, Artificial
Intelligence Graduate School Program (Sungkyunkwan University)). Kyunghyun Cho
was supported by the NYU Center for Data Science National Science Foundation
(Award 1922658) and Samsung Advanced Institute of Technology (Next Generation
Deep Learning: from pattern recognition to AI).

## References

  * Adam et al. (2018) Kalthoum Adam, Asim Baig, Somaya Al-Maadeed, Ahmed Bouridane, and Sherine El-Menshawy. 2018.  Kertas: dataset for automatic dating of ancient arabic manuscripts.  _International journal on Document Analysis and Recognition (IJDAR)_.
  * Assael et al. (2022) Yannis Assael, Thea Sommerschield, Brendan Shillingford, Mahyar Bordbar, John Pavlopoulos, Marita Chatzipanagiotou, Ion Androutsopoulos, Jonathan Prag, and Nando de Freitas. 2022.  Restoring and attributing ancient texts using deep neural networks.  _Nature_ , 603(7900):280–283.
  * Bak and Oh (2015) JinYeong Bak and Alice Oh. 2015.  Five centuries of monarchy in Korea: Mining the text of the annals of the Joseon dynasty.  In _Proceedings of the 9th SIGHUM Workshop on Language Technology for Cultural Heritage, Social Sciences, and Humanities (LaTeCH)_.
  * Bamman and Burns (2020) David Bamman and Patrick J. Burns. 2020.  Latin bert: A contextual language model for classical philology.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the NAACL_.
  * Godbole and Sarawagi (2004) Shantanu Godbole and Sunita Sarawagi. 2004.  Discriminative methods for multi-labeled classification.  In _Advances in Knowledge Discovery and Data Mining_.
  * Hossain et al. (2020) Anika Samiha Hossain, Nazia Akter, and Md. Saiful Islam. 2020.  A stylometric approach for author attribution system using neural network and machine learning classifiers.  In _Proceedings of the International Conference on Computing Advancements_.
  * Jin et al. (2020) KyoHoon Jin, JeongA Wi, KyeongPil Kang, and YoungBin Kim. 2020.  Korean historical documents analysis with improved dynamic word embedding.  _Applied Sciences_.
  * Kang et al. (2021) Kyeongpil Kang, Kyohoon Jin, Soyoung Yang, Soojin Jang, Jaegul Choo, and Youngbin Kim. 2021.  Restoring and mining the records of the Joseon dynasty via neural language modeling and machine translation.  In _Proceedings of the NAACL: Human Language Technologies_.
  * Nogueira and Cho (2019) Rodrigo Nogueira and Kyunghyun Cho. 2019.  Passage re-ranking with BERT.  _CoRR_.
  * Ouamour and Sayoud (2012) Siham Ouamour and Halim Sayoud. 2012.  authorship attribution of ancient texts written by ten arabic travelers using a smo-svm classifier.  In _International Conference on Communications and Information Technology (ICCIT)_.
  * Palladino et al. (2020) Chiara Palladino, Farimah Karimi, and Brigitte Mathiak. 2020.  Ner on ancient greek with minimal annotation.  In _https://dh2020. adho. org/_.
  * Park et al. (2020) Chanjun Park, Chanhee Lee, Yeongwook Yang, and Heuiseok Lim. 2020.  Ancient korean neural machine translation.  _IEEE Access_.
  * Reisi and Mahboob Farimani (2020) Ehsan Reisi and Hassan Mahboob Farimani. 2020.  Authorship attribution in historical and literary texts by a deep learning classifier.  _journal of Applied Intelligent Systems and Information Sciences_.
  * Robertson and Zaragoza (2009) Stephen Robertson and Hugo Zaragoza. 2009.  The probabilistic relevance framework: Bm25 and beyond.  _The Probabilistic Relevance Framework (PRF)_.
  * Sayoud and Ouamour (2017) Halim Sayoud and Siham Ouamour. 2017.  Score fusion based authorship attribution of ancient arabic texts.  In _Florida Artificial Intelligence Research Society Conference_.
  * Tian et al. (2021) Huishuang Tian, Kexin Yang, Dayiheng Liu, and Jiancheng Lv. 2021.  Anchibert: A pre-trained model for ancient chinese language understanding and generation.  In _2021 International Joint Conference on Neural Networks (IJCNN)_.
  * Vale de Gato (2015) Margarida Vale de Gato. 2015.  The collaborative anthology in the literary translation course.  _The Interpreter and Translator Trainer_.
  * Won et al. (2018) Miguel Won, Patricia Murrieta-Flores, and Bruno Martins. 2018.  Ensemble named entity recognition (ner): Evaluating ner tools in the identification of place names in historical corpora.  _Frontiers in Digital Humanities_.
  * Youden (1950) William J Youden. 1950.  Index for rating diagnostic tests.  _Cancer_ , 3(1):32–35.

## Appendix

##  Appendix A Model

Table 8 shows hyperparameter settings of our models. We used Intel(R) Xeon(R)
Silver 4114 (40 CPUs) and GeForce RTX 2080 Ti 10GB (4 GPUs) for all
experiments including training, fine-tuning, and inference.

Hyperparameter | Value
---|---
Batch Size | 32
Early Stopping Patience | 3
Hidden Size | 768
Learning Rate | 2e-5
Learning Rate Scheduler | Linear
Max Sequence Length | 512
Number of Hidden Layers | 12
Optimizer | AdamW
Vocab Size | 11270
Table 8: Model configuration

##  Appendix B HUE Dataset

###  B.1 Dataset Size

| Train | Dev | Test
---|---|---|---
CA | 330,469 | 41,309 | 41,309
TC | 330,424 | 41,303 | 41,304
NER | 385,915 | 13,417 | 13,418
SR | 169,840 | 21,570 | 21,296
Table 9: Data split in HUE dataset

###  B.2 Source Corpora

####  B.2.1 Data Collection Process

We crawl AJD 777http://sillok.history.go.kr/main/main.do, DRS
888http://sjw.history.go.kr/main.do, and DRRI
999http://kyudb.snu.ac.kr/series/main.do?item_cd=ILS from the comprehensive
database for Korean classics which are publicly available published by IKTC.
All source corpora are fully tagged with the written ages and named entities,
while their entity maps differ to each other. AJD also provides topics which
is tagged by the experts in the translation process.

####  B.2.2 Dataset Preprocessing

|  Good Examples |  Bad Examples
---|---|---
Date |  King Jeongjo 17 (1793) Feb 06 |  King Yeongjo 50 (1774) May 15
Gang |  遞承旨 徐榮輔以沈晉賢代之前望也 |  行召對于尊賢閣。 兼弼善 洪景顏。說書 李啇駿。翊贊 李應重 。
Mok |  榮輔不爲仕進政院請牌招敎以許遞前望單子入之待下批牌招察任 |  講續綱目
Gang (En) |  Royal secretary Yeong-bo Seo was replaced and Jin-hyeon Shim was appointed. |  Crown Prince held a royal lecture at the Office of Crown Prince. Kyung-ahn Hong, a lecturer for the Crown Prince, Sangjoon Lee, the second tutor of the Office of Lectures for the Crown Prince, and Eungjoon Lee, a Guard of Crown Prince attended.
Mok (En) |  When Yeong-bo Seo did not resign, the king summoned his servants through his royal secretary and ordered him to do so, saying, “wait for appointing the royal secretary among candidates and let him check his job”. |  They delivered the lecture on 《Sokgangmok (Comprehensive Mirror for Aid in Government)》.
Table 10: Good and bad examples in DRRI

Table 10 shows good and bad example in DRRI to use as summary retrieval
dataset. Bad examples mostly written from the 21 King Yeongjo to early in the
22 King Jungjo describe daily lifes of the crown prince who is King Jeongjo.
The bad example in Table 10 is depicting his study. These examples tend to
present extremely short moks which cannot be treated as summary and content,
while the offical records on administrative has much longer moks.

##  Appendix C Discussion

###  C.1 Trigram Changes Over Time

Figure 5: Trigram changes over time

Figure 5 shows the changes of trigrams over all kings. It clearly delivers the
changes of trigrams as time goes by. We can observe the same trend in either
unigram or bigram.

###  C.2 Top-5 Named Entities by Kings

| King
---
(Reigning period)
| Munjong (5 King)
---
(1450-1452)
| Seonjo (14 King)
---
(1567-1608)
| Sunjo (23 King)
---
(1790-1834)
Person | 世宗 | 2.42% | 柳成龍 | 1.08% | 南公轍 | 1.21%
金宗瑞 | 2.26% | 李德馨 | 0.81% | 金載瓚 | 1.09%
李季甸 | 1.97% | 尹斗壽 | 0.72% | 李時秀 | 0.85%
皇甫仁 | 1.82% | 李恒福 | 0.65% | 沈象奎 | 0.83%
鄭苯 | 1.78% | 李元翼 | 0.57% | 李相璜 | 0.77%
Person (En) | King Sejong | (1392-1397) | Seongryong Ryu | (1542-1607) | Gongcheol Nam | (1760-1840)
Jongseo Kim | (1383-1453) | Deokhyeong Lee | (1561-1613) | Jaechan Kim | (1746-1827)
Kyejeon Lee | (1404-1459) | Dushu Yun | (1533-1601) | Sisu Lee | (1745-1821)
Boin Hwang | (1387-1453) | Hangbok Lee | (1556-1618) | Sangkyu Shim | (1766-1838)
Bun Jeong | (1394-1454) | Weonik Lee | (1547-1634) | Sanghwang Lee | (1763-1841)
Location | 平安 | 4.40% | 平壤 | 2.10% | 春塘臺 | 2.02%
咸吉 | 3.24% | 朝鮮 | 1.93% | 漢城府 | 1.19%
黃海 | 2.28% | 京畿 | 1.17% | 仁政殿 | 1.71%
輝德殿 | 2.21% | 全羅 | 1.76% | 平安 | 1.51%
京畿 | 2.10% | 慶尙 | 1.64% | 景慕宮 | 1.49%
Location (En) | Pyong-an | (Province) |  Pyongyang |  (Province) | Chundangdae | (Site)
Hamgyong | (Province) | Joseon | (Country) | Hanseong Magistracy | (Province)
Hwanghae |  (Province) | Gyeonggi | (Province) | Injeongjeon | (Hall)
Hwideokjeon | (Hall) |  Jeolla |  (Province) | Pyong-an | (Province)
Gyeonggi | (Province) |  Gyeongsang |  (Province) | Gyeongmogung | (Palace)

Table 11: Top-5 named entities in 5, 14, and 23 kings

Table 11 shows top-5 person and location named entities in three King reigns.
All person entities except King Sejong, the most frequent entity in Munjong,
are officials in the reign periods, which is different by time changes. In
contary, all location entities are the name of place, palace or site, showing
some entities overlap among kings.

##  Appendix D Experimental Results

| Acc | F1 | Pre | Rec | QWK
---|---|---|---|---|---
BERT not pretrained | 56.59 | 54.26 | 55.08 | 56.59 | 87.09
AnchiBERT Tian et al. (2021) | 76.31 | 75.74 | 76.45 | 76.31 | 93.99
mBERT Devlin et al. (2019) | 76.02 | 75.29 | 75.80 | 76.02 | 93.76
AnchiBERT + AJD/DRS | 79.50 | 79.33 | 80.06 | 79.50 | 95.46
mBERT + AJD/DRS | 77.99 | 77.78 | 78.92 | 77.99 | 95.04
Table 12: Evaluation results for our PLMs on chronological attribution

For CA, we measure the Quadratic Weighted Kappa score (QWK score) as metrics
that treat each king label hierarchically.

| Major (4 classes) | Minor (106 classes)
---|---|---
| Acc | F1 | Pre | Rec | Ham | Acc | F1 | Pre | Rec | Ham
BERT not pretrained | 68.52 | 68.91 | 74.82 | 70.65 | 79.04 | 26.48 | 61.29 | 54.86 | 84.59 | 42.45
AnchiBERT Tian et al. (2021) | 68.99 | 85.81 | 85.24 | 88.44 | 83.00 | 31.47 | 69.30 | 61.63 | 90.26 | 51.63
mBERT Devlin et al. (2019) | 56.78 | 79.59 | 77.60 | 87.78 | 73.64 | 32.85 | 70.01 | 62.17 | 90.84 | 54.28
AnchiBERT + AJD/DRS | 70.48 | 88.33 | 86.61 | 92.71 | 84.80 | 31.77 | 72.07 | 64.87 | 91.24 | 52.25
mBERT + AJD/DRS | 69.15 | 87.13 | 85.48 | 91.36 | 83.81 | 33.96 | 71.28 | 63.50 | 91.56 | 55.91

Table 13: Evaluation results for our PLMs on topic classification

Since TC is a multi-label classification task whose example might have
multiple labels as the answer, we measure Hamming score along with accuracy.
In this case, accuracy is the exact match score, and the hamming score is the
accuracy of subset matched, $\mid T\cap P\mid/\mid T\cup P\mid$, where $T$ is
set of true labels and $P$ is set of predicted labels Godbole and Sarawagi
(2004). For the evaluation results in Table 13, we find and set the best
threshold to each label by Youden’s index. All pretrained models outperform
BERT without pretraining, and two LMs re-trained on hanja documents show the
best performances.

| Overall | Person | Location
---|---|---|---
| Acc | F1 | Pre | Rec | F1 | Pre | Rec | F1 | Pre | Rec
BERT not pretrained | 98.67 | 89.40 | 90.58 | 88.25 | 92.13 | 93.98 | 90.36 | 87.10 | 87.76 | 86.46
AnchiBERT | 98.72 | 90.30 | 90.98 | 89.62 | 93.13 | 94.47 | 91.83 | 87.91 | 88.08 | 87.75
mBERT | 98.52 | 88.57 | 89.52 | 87.64 | 91.63 | 93.60 | 89.74 | 86.02 | 86.18 | 85.85
AnchiBERT + AJD/DRS | 98.76 | 90.42 | 91.31 | 89.55 | 93.28 | 94.77 | 91.84 | 88.01 | 88.43 | 87.59
mBERT + AJD/DRS | 98.69 | 90.16 | 90.36 | 89.95 | 92.83 | 93.51 | 92.17 | 87.90 | 87.74 | 88.06

Table 14: Evaluation results for our PLMs on named entity recognition | MRR | Top-1 | Top-10
---|---|---|---
BM25 | 29.87 | 25.58 | 33.98
BERT not pretrained | 52.85 | 99.20 | 99.64
AnchiBERT Tian et al. (2021) | 67.92 | 99.20 | 99.85
mBERT Devlin et al. (2019) | 67.06 | 99.32 | 99.50
AnchiBERT + AJD/DRS | 74.29 | 99.64 | 99.91
mBERT + AJD/DRS | 73.88 | 99.44 | 99.59
Table 15: Evaluation results for our PLMs on summary retrieval


The amplitude $|A(\tau^{\prime})\\}$, which depends on the optical depth
$\tau^{\prime}$ of the point source, can be determined from the boundary
condition that no thermal radiation generated by the cloud comes into the top
or bottom.

The value of the Green’s function (276) for observation points at the top of
the cloud, where $\tau=\tau_{c}$, is

$\displaystyle|G(\tau_{c},\tau^{\prime})\\}$ $\displaystyle=$
$\displaystyle|G^{\\{\infty\\}}(\tau_{c}-\tau^{\prime})\\}+|\Delta
G(\tau_{c},\tau^{\prime})\\}$ (282) $\displaystyle=$
$\displaystyle\hat{\kappa}_{\bf u}e^{-\hat{\kappa}_{\bf
u}(\tau_{c}-\tau^{\prime})}|0)+\left[\mathcal{L}_{\bf d}+e^{-\hat{\kappa}_{\bf
u}\tau_{c}}\right]|A(\tau^{\prime})\\}.$

For observation points at the bottom of the cloud, where $\tau=0$, we find

$\displaystyle|G(0,\tau^{\prime})\\}$ $\displaystyle=$
$\displaystyle|G^{\\{\infty\\}}(-\tau^{\prime})\\}+|\Delta
G(0,\tau^{\prime})\\}$ (283) $\displaystyle=$ $\displaystyle-\hat{\kappa}_{\bf
d}e^{\hat{\kappa}_{\bf d}\tau^{\prime}}|0)+\left[e^{\hat{\kappa}_{\bf
d}\tau_{c}}+\mathcal{L}_{\bf u}\right]|A(\tau^{\prime})\\}.$

Using (282) with the boundary condition, $\mathcal{M}_{\bf
d}|G(\tau_{c},\tau^{\prime})\\}=\hat{0}$, of no downward internally generated
thermal radiation at the top of the cloud, and using (283) with the boundary
condition, $\mathcal{M}_{\bf u}|G(0,\tau^{\prime})\\}=\hat{0}$ of no
internally generated upward thermal radiation at the bottom of the cloud, we
find

$\displaystyle\hat{0}$ $\displaystyle=$ $\displaystyle\mathcal{C}_{\bf
du}\hat{\kappa}_{\bf u}e^{-\hat{\kappa}_{\bf
u}(\tau_{c}-\tau^{\prime})}|0)+\left[\mathcal{C}_{\bf dd}+\mathcal{C}_{\bf
du}e^{-\hat{\kappa}_{\bf u}\tau_{c}}\right]|A(\tau^{\prime})\\},$ (284)
$\displaystyle\hat{0}$ $\displaystyle=$ $\displaystyle-\mathcal{C}_{\bf
ud}\hat{\kappa}_{\bf d}e^{\hat{\kappa}_{\bf
d}\tau^{\prime}}|0)+\left[\mathcal{C}_{\bf ud}e^{\hat{\kappa}_{\bf
d}\tau_{c}}+\mathcal{C}_{\bf uu}\right]|A(\tau^{\prime})\\}.$ (285)

Summing the left and right sides of equations (284) and (285) we find

$\displaystyle\hat{0}$ $\displaystyle=$ $\displaystyle\left[\mathcal{C}_{\bf
du}\hat{\kappa}_{\bf u}e^{-\hat{\kappa}_{\bf
u}(\tau_{c}-\tau^{\prime})}-\mathcal{C}_{\bf ud}\hat{\kappa}_{\bf
d}e^{\hat{\kappa}_{\bf d}\tau^{\prime}}\right]|0)$ (286) $\displaystyle+$
$\displaystyle\left[\mathcal{C}_{\bf dd}+\mathcal{C}_{\bf
du}e^{-\hat{\kappa}_{\bf u}\tau_{c}}+\mathcal{C}_{\bf ud}e^{\hat{\kappa}_{\bf
d}\tau_{c}}+\mathcal{C}_{\bf uu}\right]|A(\tau^{\prime})\\}.$

Introducing the matrix

$\mathcal{R}(\tau^{\prime})=-\mathcal{C}_{\bf du}e^{-\hat{\kappa}_{\bf
u}(\tau_{c}-\tau^{\prime})}+\mathcal{C}_{\bf ud}e^{\hat{\kappa}_{\bf
d}\tau^{\prime}},$ (287)

and noting the definition (206) of $\mathcal{I}$, we rewrite (286) as

$\frac{\partial}{\partial\tau^{\prime}}\mathcal{R}(\tau^{\prime})|0)=\mathcal{I}|A(\tau^{\prime})\\}.$
(288)

Since $\mathcal{I}$ is independent of $\tau^{\prime}$ we can solve (288) for
$|A(\tau^{\prime})\\}$ to find

$|A(\tau^{\prime})\\}=\frac{\partial}{\partial\tau^{\prime}}\mathcal{I}^{-1}\mathcal{R}(\tau^{\prime})|0)$
(289)

Substituting (289) into (281) we find

$|\Delta
G(\tau,\tau^{\prime})\\}=\frac{\partial}{\partial\tau^{\prime}}\mathcal{U}(\tau)\mathcal{I}^{-1}\mathcal{R}(\tau^{\prime})|0).$
(290)

Adding (277) to (290) gives an explicit expression for
$|G(\tau,\tau^{\prime})\\}$, the total Green’s function (276) of the cloud.

$\displaystyle|G(\tau,\tau^{\prime})\\}=\delta(\tau-\tau^{\prime})|0)$
$\displaystyle+\frac{\partial}{\partial\tau^{\prime}}\bigg{[}H(\tau-\tau^{\prime})e^{-\hat{\kappa}_{\bf
u}(\tau-\tau^{\prime})}-H(\tau^{\prime}-\tau)e^{-\hat{\kappa}_{\bf
d}(\tau-\tau^{\prime})}+\mathcal{U}(\tau)\mathcal{I}^{-1}\mathcal{R}(\tau^{\prime})\bigg{]}|0).$
(291)

Figure 19: Thermal emission, computed with (294) for a forward-scattering
cloud of optical depth $\tau_{c}=8$. The cloud is isothermal with a constant
value of the Planck intensity $B=1$. The radiation is modeled with $2n=10$
streams. The lengths of the blue rays are proportional to the intensities of
the thermal radiation,
$T(\mu_{i},\tau)=w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|T(\tau)\\}$, generated in
the cloud. The forward-scattering phase function is constructed from the first
$2p=6$ Legendre polynomials, like that of Fig. 4. The single-scattering albedo
is $\tilde{\omega}=0.9$. At the center of the cloud, at an optical depth
$\Delta\tau=4$ from both the top and bottom, the thermal radiation is nearly
isotropic, but slightly more intense in horizontal than vertical directions.
At the top and bottom of the cloud, marked by the horizontal black lines,
there is only outgoing, “limb-darkened” thermal radiation.

Figure 20: Like Fig. 19, but for an isothermal cloud with optical depth
$\tau_{c}=8$, located above a warmer black surface with Planck intensity
$1.2B$ that is 20% larger than the Planck intensity $B$ of the cloud. The
intensity in the cloud is a combination of “blue” thermal emission by cloud
particulates, computed with (294), and attenuated “red” thermal radiation from
the warmer surface, calculated with (210). The total intensity is the sum of
the lengths of the two colored segments. The upward radiation from the surface
has the incident intensity $|I^{\\{\rm in\\}}\\}=\mathcal{M}_{\bf u}|0)1.2B$
and is Lambertian. At increasing optical depths above the surface, the slant
rays of the red surface radiation are more attenuated than the more vertical
radiation, which leads to increasing limb darkening with altitude. Very little
“red” surface radiation emerges from the cloud top, where the outgoing
radiation is almost all “blue” thermal radiation generated in the cloud. A
relatively small amount of “red” surface radiation is diffusely reflected from
the bottom of the cloud.

###  9.1 Isothermal emitting cloud

If the cloud is isothermal, the Planck intensity $B(\tau^{\prime})$ can be
taken to have the $\tau^{\prime}$-independent value $B$. Then we can use (273)
with (291) to write the thermal intensity in the cloud as

$\displaystyle|T(\tau)\\}=B\int_{0}^{\tau_{c}}d\tau^{\prime}|G(\tau,\tau^{\prime})\\}$
(292) $\displaystyle=$
$\displaystyle|B\\}+\bigg{[}H(\tau-\tau^{\prime})e^{-\hat{\kappa}_{\bf
u}(\tau-\tau^{\prime})}-H(\tau^{\prime}-\tau)e^{-\hat{\kappa}_{\bf
d}(\tau-\tau^{\prime})}+\mathcal{U}(\tau)\mathcal{I}^{-1}\mathcal{R}(\tau^{\prime})\bigg{]}_{\tau^{\prime}=0}^{\tau^{\prime}=\tau_{c}}|B\\}$
$\displaystyle=$ $\displaystyle\bigg{[}\hat{1}-e^{-\hat{\kappa}_{\bf
d}(\tau-\tau_{c})}+\mathcal{U}(\tau)\mathcal{I}^{-1}\mathcal{R}(\tau_{c})-e^{-\hat{\kappa}_{\bf
u}\tau}-\mathcal{U}(\tau)\mathcal{I}^{-1}\mathcal{R}(0)\bigg{]}|B\\}$
$\displaystyle=$
$\displaystyle\bigg{[}\hat{1}-\mathcal{U}(\tau)+\mathcal{U}(\tau)\mathcal{I}^{-1}\\{\mathcal{R}(\tau_{c})-\mathcal{R}(0)\\}\bigg{]}|B\\}.$

From (287), (206) and (170) we find

$\displaystyle\mathcal{R}(\tau_{c})-\mathcal{R}(0)$ $\displaystyle=$
$\displaystyle-\mathcal{C}_{\bf du}+\mathcal{C}_{\bf ud}e^{\hat{\kappa}_{\bf
d}\tau_{c}}+\mathcal{C}_{\bf du}e^{-\hat{\kappa}_{\bf
u}\tau_{c}}-\mathcal{C}_{\bf ud}$ (293) $\displaystyle=$
$\displaystyle\mathcal{I}-\hat{1}.$

Substituting (293) into (292) we find that the thermal radiation generated by
the cloud is

$|T(\tau)\\}=\bigg{[}\hat{1}-\mathcal{U}(\tau)\mathcal{I}^{-1}\bigg{]}|B\\}.$
(294)

From (294) we see that the thermal outgoing radiation from the cloud is

$\displaystyle|T^{\\{\rm out\\}}\\}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{\bf d}|T(0)\\}+\mathcal{M}_{\bf u}|T(\tau_{c})\\}$
(295) $\displaystyle=$ $\displaystyle\left[\mathcal{M}_{\bf
d}+\mathcal{M}_{\bf u}-\\{\mathcal{M}_{\bf d}\,\mathcal{U}(0)+\mathcal{M}_{\bf
u}\,\mathcal{U}(\tau_{c})\\}\mathcal{I}^{-1}\right]|B\\}$ $\displaystyle=$
$\displaystyle\left[\hat{1}-\mathcal{O}\mathcal{I}^{-1}\right]|B\\}$
$\displaystyle=$ $\displaystyle\mathcal{E}|B\\}.$

The outgoing matrix $\mathcal{O}$ was defined by (208). According to (213),
the scattering matrix of the cloud is
$\mathcal{S}=\mathcal{O}\mathcal{I}^{-1}$. So from (297) we see that the
emissivity matrix is

$\mathcal{E}=\hat{1}-\mathcal{O}\mathcal{I}^{-1}=\hat{1}-\mathcal{S}.$ (296)

Eq. (296) is Kirchhoff’s law (185). As a consistency check, we use (294) and
(206) to show that the thermal radiation coming into the cloud is

$\displaystyle|T^{\\{\rm in\\}}\\}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{\bf u}|T(0)\\}+\mathcal{M}_{\bf d}|T(\tau_{c})\\}$
(297) $\displaystyle=$ $\displaystyle\left[\mathcal{M}_{\bf
u}+\mathcal{M}_{\bf d}-\\{\mathcal{M}_{\bf
u}\mathcal{\,}{U}(0)+\mathcal{M}_{\bf
d}\,\mathcal{U}(\tau_{c})\\}\mathcal{I}^{-1}\right]|B\\}$ $\displaystyle=$
$\displaystyle\left[\hat{1}-\mathcal{I}\mathcal{I}^{-1}\right]|B\\}$
$\displaystyle=$ $\displaystyle\hat{0}.$

##  10 Summary

Here we summarize the most important differences between the results of our
$2n$-stream model and previous work. Additional interesting and useful results
will be discussed in subsequent papers.

We consider axially symmetric variables, like the intensity $I(\mu,\tau)$,
that depend on the direction cosine $\mu$ of the radiation, and on the optical
depth $\tau$ above the bottom of a cloud. Using Dirac notation, similar to
that used for the analysis of Schrödinger’s equation in quantum mechanics, we
write the equation of transfer in the formally simple, but completely general
form (69). The propagation of radiation is controlled by the $2n\times 2n$
exponentiation-rate matrix $\hat{\kappa}$ of (70) that is analogous to the
quantum mechanical Hamiltonian operator $H$. The intensity is described by an
abstract state vector $|I(\tau)\\}$ that depends on the optical depth $\tau$
in a similar way to the dependence of a spin wave function $|\psi(t)\rangle$
on the time $t$.

As exemplified in (128), we represent radiative variables as a superposition
of $2n$ basis vectors, the right stream bases $|\mu_{i})$ and the left stream
bases $(\kern-2.59995pt(\mu_{i}|$ of (102) and (103). The stream indices are
$i=1,2,3,\ldots,2n$. The streams have nominal zenith angles
$\theta_{i}=\cos^{-1}\mu_{i}$. According to (73), the Gauss-Legendre direction
cosine $\mu_{i}$ is a root of the Legendre polynomial $P_{2n}$, that is,
$P_{2n}(\mu_{i})=0$. Representative projections $\langle\mu|\mu_{i})$ of the
stream bases onto continuous direction-cosine states $\langle\mu|$ of (26) are
shown in Fig. 1. The maximum values of $|\langle\mu|\mu_{i})|$ occur for
$\mu\approx\mu_{i}$. So the stream bases $|\mu_{i})$ represent radiation
propagating with directions close to that of the zenith angle
$\theta_{i}=\cos^{-1}\mu_{i}$.

As in quantum mechanics, we think of bases as eigenvectors of “observables.”
In (112) we show that the stream bases $|\mu_{i})$ are right eigenvectors of
the $2n\times 2n$ direction-cosine matrix $\hat{\mu}$ of (36), or of its
inverse, the direction-secant matrix $\hat{\varsigma}$ of (51). The left
eigenvectors of $\hat{\mu}$ are denoted by $(\kern-2.59995pt(\mu_{i}|$. The
double left parenthesis is a reminder that the left eigenvectors are not
Hermitian conjugates of the right eigenvectors. Many radiative-transfer
variables in $2n$ space can be conveniently represented with non-Hermitian
matrices. So it is not always possible to write left and right eigenvectors as
Hermitian conjugate pairs.

Radiative transfer is parameterized by the single scattering-albedo
$\tilde{\omega}$ and by the scattering phase function $p(\mu,\mu^{\prime})$.
The fraction of photons scattered by a collision is $\tilde{\omega}$. The
fraction absorbed and converted to heat is $1-\tilde{\omega}$. The fraction of
photons of initial direction cosine $\mu^{\prime}$ that is scattered into
direction cosines between $\mu$ and $\mu+d\mu$ is
$p(\mu,\mu^{\prime})\,d\mu/2$.

We have introduced a family of phase functions,
$p(\mu,\mu^{\prime}=1)=\varpi^{\\{p\\}}(\mu)$ of (141) which are
superpositions of the first $2p$ Legendre polynomials, and which produce the
maximum possible forward scattering. For a $2n$-stream model one must have
$p\leq n$. The maximum-forward-scattering phase function
$\varpi^{\\{p\\}}(\mu)$ satisfies the constraints (5) and (6). The forward-
scattering values are $\varpi^{\\{p\\}}(\mu=1)=p(p+1)$. A representative
forward-scattering phase function $\varpi^{\\{3\\}}(\mu,\mu^{\prime})$ is
shown in Fig. 4, and others are shown in Fig. 6.

The homogeneous equation of transfer (198) has $2n$ independent solutions,
$|\lambda_{i})e^{-\kappa_{i}\tau}$. The penetration modes $|\lambda_{i})$,
defined by (151), are eigenfunctions of the exponentiation-rate matrix
$\hat{\kappa}$, or of its inverse, $\hat{\lambda}=\hat{\kappa}^{-1}$, the
penetration-length matrix. Representative projections
$\langle\mu|\lambda_{i})$ of the penetration modes onto the direction-cosine
states $\langle\mu|$ are shown in Fig. 8 for a weakly absorbing cloud with
$1-\tilde{\omega}=0.001$. For negligible scattering, $\tilde{\omega}\to 0$, we
noted in (171) that the penetration modes are identical to the stream bases,
$|\lambda_{i})=|\mu_{i})$, and the penetration lengths are the same as the
Gauss-Legendre cosines, $\lambda_{i}=\mu_{i}$. For any value of
$\tilde{\omega}<0$, the directional penetration modes $|\lambda_{i})$, with
$i=2,3,\dots,2n-1$, are qualitatively similar to the corresponding stream
bases $|\mu_{i})$. But for weak absorption, the upward quasi-isotropic mode
$|\lambda_{2n})$, shown in the left panel of Fig. 8, has very little
dependence on $\mu$, and it differs drastically from the stream basis
$|\mu_{2n})$ on the left panel of Fig. 1, which is strongly peaked near
$\mu=1$. The downward quasi-isotropic mode $|\lambda_{1})$ differs from
$|\mu_{1})$ in an analogous way. The quasi-isotropic modes represent the
familiar, nearly uniform brightness seen from all directions on a foggy day or
from an aircraft flying through thick, sunlit clouds.

Representative penetration lengths, $\lambda_{i}=1/\kappa_{i}$, the inverses
of the exponentiation rates $\kappa_{i}$, are plotted in Fig. 7 as a function
of the single-scattering albedo $\tilde{\omega}$. As $\tilde{\omega}\to 1$,
the penetration lengths $\lambda_{i}$ separate into two groups. For the
directional modes, with $i=2,3,\ldots,2n-1$, the lengths $|\lambda_{i}|$ are
only slightly larger than the corresponding Gauss-Legendre cosines.
$|\mu_{i}|$. The quasi-isotropic mode $|\lambda_{2n})$, and its mirror image,
$|\lambda_{1})$ have divergent penetration lengths, as one can see by
inspection of Fig. 7. According to (174), as $\tilde{\omega}\to 1$ the quasi-
isotropic penetration lengths are very nearly given by
$\lambda_{2n}=-\lambda_{1}\to\big{[}3(1-\tilde{\omega}p_{1})(1-\tilde{\omega})\big{]}^{1/2}$.

The optical-depth evolution operator $\mathcal{U}(\tau)$ of (200) is the sum
of exponentially growing or decaying terms of the form
$e^{-\kappa_{i}(\tau-\tau_{i})}|\lambda_{i})(\kern-2.59995pt(\lambda_{i}|$. As
indicated by (201), for radiative transfer it is convenient to choose
reference optical depths $\tau_{i}=0$, at the bottom of the cloud, for
decaying modes with $\kappa_{i}>0$. For growing modes, with $\kappa_{i}<0$,
the reference optical depths, $\tau_{i}=\tau_{c}$, are most conveniently
chosen at the top of the cloud. According to (210), the intensity
$|I(\tau)\\}$ at an optical depth $\tau$ inside a cloud is given in terms of
the incident intensity $|I^{\\{\rm in\\}}\\}$ at the top and bottom of the
cloud by $|I(\tau)\\}=\mathcal{U}(\tau)\mathcal{I}^{-1}|I^{\\{\rm in\\}}\\}$,
where the incoming matrix $\mathcal{I}$ was given by (206).

We represent the overall radiation transfer in non-emissive clouds by the
scattering matrix $\mathcal{S}$ of (213). According to (212), the outgoing
intensity, $|I^{\\{\rm out\\}}\\}$, from the top and bottom of the cloud is
related to the incoming intensity $|I^{\\{\rm in\\}}\\}$ by $|I^{\\{\rm
out\\}}\\}=\mathcal{S}|I^{\\{\rm in\\}}\\}$. The scattering matrices
$\mathcal{S}^{(1)}$ and $\mathcal{S}^{(2)}$ of two clouds in series can be
represented by a single scattering matrix $\mathcal{S}$, given by the addition
formula of (252).

We consider not only the scattering of external radiation by clouds but also
the thermal emission of radiation by warm molecules or particulates inside
clouds. The Green’s function (or point-response function)
$G(\tau,\tau^{\prime})$ of (291) can be used to calculate the thermally
emitted intensity $|T(\tau)\\}$ within a non-isothermal cloud with (273). For
the simple case of an isothermal cloud with a constant Planck intensity
$|B\\}$, (295) shows that the thermally generated intensity going out of the
cloud is $|T^{\\{\rm out\\}}\\}=\mathcal{E}|B\\}$, where $\mathcal{E}$ is the
emissivity matrix of the cloud. In agreement with Kirchhoff’s law, (296) shows
that the sum of the isothermal emissivity matrix and the scattering matrix is
the unit matrix, $\mathcal{E}+\mathcal{S}=\hat{1}$.

Fig. 16 compares diffuse reflection and transmission for pure isotropic
scattering in clouds of various optical depths $\tau_{c}$. The top panels of
Fig. 16 reproduce Fig. 21 of Chandrasekhar[2]. The bottom panels are the
results of a 32-stream calculation with $\tilde{\omega}=1-10^{-6}$ and with an
input stream $|I^{\\{\rm in\\}}\\}=|\mu_{23})$, where $\mu_{23}=0.5877$ is the
Gauss-Legendre cosine that is closest to the input cosine $\mu=0.6$ used by
Chandrasekhar. In our Fig. 17 we show a similar comparison to Chandrasekhar’s
Fig. 22 where the single-scattering albedo was $\tilde{\omega}=0.9$. Our
results can hardly be distinguished from those of Chandrasekhar, who used “$X$
and $Y$ functions” of §61(1) and §61(2), the solutions of complicated integral
equations.

The $2n$ values of intensity,
$I(\mu_{i},\tau)=w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$, can be used
to make physically instructive plots like Fig. 10, which shows how part the
nearly vertical radiation incident on the bottom of a cloud is diffusely
reflected and part is transmitted. For optically thick clouds like that of
Fig. 10, with $\tau_{c}=8$, but with very small single-scattering absorption
$1-\tilde{\omega}=10^{-2}$, appreciable transmission (12%) is possible because
of multiple scattering. Multiple scattering nearly “isotropizes” the intensity
near the center of the cloud. At the top of the cloud, there are upward
transmitted streams, but no downward streams. At the bottom of the cloud there
is a strong incident stream $|I^{\\{\rm in\\}}\\}=|\mu_{2n})$, and relatively
weak diffusely reflected streams $I(\mu_{j},\tau=0)$ with $\mu_{j}<0$.

Fig. 11 shows what happens when the single-scattering absorption is reduced by
a factor of 100 from that of Fig. 10, to $1-\tilde{\omega}=10^{-4}$. For such
small single-scattering absorption, the absorption in the cloud is also
negligible and the flux is nearly constant from the bottom to the top of the
cloud. The formalism outlined in this paper does not work if there is no
absorption at all and $\tilde{\omega}=1$. As we will show in a subsequent
paper, minor modifications in the computational steps outlined above allow us
to calculate radiative tranfer for $\tilde{\omega}=1$, Chandrasekhar’s
“conservative scattering.”

Other interesting calculations are summarized in Fig. 13, which shows how the
transmission of radiation from the bottom to the top of a cloud with
negligibly small absorption $1-\tilde{\omega}\ll 1$ depends on the optical
depth $\tau_{c}$ of the cloud and on the nature of the scattering phase
function $p(\mu,\mu^{\prime})$. Not suprisingly, the transmission is greatest
for forward scattering and least for backward scattering.

Fig. 15, for an optically thick cloud, shows the cloud albedo,
$\langle\Omega\rangle$, of (224). The cloud albedo is the fraction of the
vertical flux reflected from the bottom of a cloud, as opposed to being
absorbed. The cloud albedo depends on the single-scattering albedo
$\tilde{\omega}$, on the nature of the scattering phase function
$p(\mu,\mu^{\prime})$, and on the angular distribution of the incident
intensity $|I^{\\{\rm in\\}}\\}$. If the cloud were not optically thick, its
albedo would also depend on the optical thickness $\tau_{c}$.

Finally, Fig. 20 shows how blackbody thermal flux emitted from a warm surface
propagates into a cooler, isothermal cloud with single-scattering albedo
$\tilde{\omega}=0.9$, a forward scattering phase function,
$\varpi^{\\{3\\}}(\mu)$, and an optical depth $\tau_{c}=8$. Most of the
surface radiation is absorbed and replaced by thermal flux emitted by the
cloud. A small fraction of the radiation is diffusely reflected back into the
surface, and a still smaller fraction is transmitted out of the top of the
cloud, where most of the emergent radiation has been thermally generated by
cloud particulates.

It is easy to use modern mathematical software like Matlab with this new
$2n$-stream formalism to calculate a variety of radiative-transfer quantities,
including diffuse reflection or transmission, albedos, thermal emission, etc.
A few dozen lines of code is usually sufficient. It is easy to go beyond
isotropic or Rayleigh-scattering phase functions of traditional radiative
transfer calculations to much more realistic and complicated ones. The family
of phase functions $\varpi^{\\{p\\}}(\mu)$ of (141) allows one to model
radiation transfer with large forward-scattering amplitudes, like the Mie
scattering of sunlight in Earth’s clouds. We believe that the $2n$-stream
formalism discussed here will be a useful new tool for radiative-transfer
calculations.

##  Appendix A Phase Functions for Maximum Forward Scattering

We turn now to the proof that the algebraic formula of (141) for
$\varpi^{\\{p\\}}(\mu)$ gives the maximum possible forward scattering for a
phase function constructed from the first $2p=2,4,6,\ldots$ Legendre
polynomials. We assume that the maximum forward-scattering phase function can
be written as

$\varpi^{\\{p\\}}(\mu)=(1+\mu)\phi^{2}(\mu).$ (298)

Here $\phi(\mu)$ is a polynomial of degree $p-1$ which we expand in Legendre
polynomials as

$\displaystyle\phi(\mu)$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{p-1}\phi_{l}\,g_{l}\,P_{l}(\mu).$ (299)

The coefficient $g_{l}$ is a diagonal element of the statistical weight matrix
$\hat{g}$, defined by (164). To facilitate further discussion we introduce a
statistical-weight column vector $|g\rangle$ and row vector $\langle g|$ with
the elements

$|g\rangle=\left[\begin{array}[]{c}g_{0}\\\ g_{1}\\\ g_{2}\\\ \vdots\\\
g_{p-1}\end{array}\right],\quad\hbox{and}\quad\langle
g|=[g_{0},\,g_{1},\,g_{2},\,\cdots,g_{p-1}].$ (300)

Using (299) we write the forward-scattering phase (298) as

$\displaystyle\varpi^{\\{p\\}}(1)$ $\displaystyle=$ $\displaystyle
2\phi^{2}(1)=2\left(\sum_{l=0}^{p-1}\phi_{l}g_{l}\right)^{2}$ (301)
$\displaystyle=$ $\displaystyle
2\sum_{ll^{\prime}}\phi_{l}q_{ll^{\prime}}\phi_{l^{\prime}}.$

In (301) we have introduced the $p\times p$ symmetric matrix $\hat{q}$,
defined by

$\hat{q}=|g\rangle\langle
g|,\quad\hbox{or}\quad(\kern-2.59995pt(l|\hat{q}|l^{\prime})=q_{ll^{\prime}}=g_{l}g_{l^{\prime}}.$
(302)

As shown in (6), the area of the phase function $\varpi^{\\{p\\}}(\mu)$ is
constrained to be

$\int_{-1}^{1}d\mu\,\varpi^{\\{p\\}}(\mu)=\int_{-1}^{1}d\mu\,(1+\mu)\sum_{ll^{\prime}}g_{l}\phi_{l}P_{l}(\mu)g_{l^{\prime}}\phi_{l^{\prime}}P_{l^{\prime}}(\mu)=2,$
(303)

or

$1=\sum_{ll^{\prime}}\phi_{l}h_{ll^{\prime}}\phi_{l^{\prime}}.$ (304)

Here the elements $h_{ll^{\prime}}$ of the constraint matrix are

$h_{ll^{\prime}}=g_{l}\delta_{ll^{\prime}}+g_{l}\hat{\mu}_{ll^{\prime}},\quad\hbox{or}\quad\hat{h}=\hat{g}+\hat{g}\hat{\mu}=\hat{g}(\hat{1}+\hat{\mu}).$
(305)

The elements of the direction-cosine matrix $\hat{\mu}$ were given by (17) or
(36), and they can be used to write elements of the constraint matrix
$\hat{h}$ of (305) as

$h_{ll^{\prime}}=(\kern-2.59995pt(l|\hat{h}|l^{\prime})=\left[\begin{array}[]{lllllll}1&1&0&0&0&0&\cdots\\\
1&3&2&0&0&0&\cdots\\\ 0&2&5&3&0&0&\cdots\\\ 0&0&3&7&4&0&\cdots\\\
0&0&0&4&9&5&\cdots\\\ 0&0&0&0&5&11&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\ \end{array}\right]$ (306)

We use a Lagrangian multiplier $m$ to find the coefficients $\phi_{l}$ that
maximize the forward-scattering phase $\varpi^{\\{p\\}}(1)$ of (301), subject
to the constraint (304). Then we can vary the amplitudes $\phi_{k}$
independently to find local maxima, minima or saddle points of the function

$M(\phi_{0},\phi_{1},\ldots,\phi_{n-1})=\varpi^{\\{p\\}}(1)-m\sum_{ll^{\prime}}\phi_{l}h_{ll^{\prime}}\phi_{l^{\prime}}.$
(307)

For $\phi_{k}$’s at a local maximum of $M$, we must have

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\phi_{k}}M$ (308) $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\phi_{k}}\sum_{ll^{\prime}}\phi_{l}\left(2q_{ll^{\prime}}-mh_{ll^{\prime}}\right)\phi_{l^{\prime}}$
$\displaystyle=$ $\displaystyle
2\sum_{l^{\prime}}\left(2q_{kl^{\prime}}-mh_{kl^{\prime}}\right)\phi_{l^{\prime}},\quad\hbox{for}\quad
k=0,1,2,\ldots,2p-1.$

In writing the last line of (308) we made use of the symmetries
$q_{ll^{\prime}}=q_{l^{\prime}l}$ and $h_{ll^{\prime}}=h_{l^{\prime}l}$. We
can write the last line of (308) as

$2\sum_{l^{\prime}}q_{kl^{\prime}}\phi_{l^{\prime}}=m\sum_{l^{\prime}}h_{kl^{\prime}}\phi_{l^{\prime}},\quad\hbox{or}\quad
2|g\rangle\langle g|\phi)=m\hat{h}|\phi).$ (309)

Here we have introduced the column vector of multipole amplitudes

$|\phi)=\left[\begin{array}[]{c}\phi_{0}\\\ \phi_{1}\\\ \phi_{2}\\\ \vdots\\\
\phi_{p-1}\end{array}\right].$ (310)

We can multiply the left version of (309) by $\phi_{k}$, sum over $k$, and use
(301) with (304) to find

$\displaystyle\varpi^{\\{p\\}}(1)$ $\displaystyle=$ $\displaystyle
2\sum_{kl^{\prime}}\phi_{k}q_{kl^{\prime}}\phi_{l^{\prime}}$ (311)
$\displaystyle=$ $\displaystyle
m\sum_{kl^{\prime}}\phi_{k}h_{kl^{\prime}}\phi_{l^{\prime}}$ $\displaystyle=$
$\displaystyle m.$

The maximum value $\varpi^{\\{p\\}}(1)$ of the phase function is the same as
the value of $m$, the Lagrange multiplier. Multiplying both sides of the
second, matrix equation of (309) by $\hat{h}^{-1}$ we find

$\hat{m}|\phi\rangle=m|\phi\rangle,\quad\hbox{where}\quad\hat{m}=|f\rangle\langle
g|.$ (312)

The column vector $|f\rangle$ of (312) is

$\displaystyle|f\rangle$ $\displaystyle=$ $\displaystyle
2\hat{h}^{-1}|g\rangle.$ (313)

To facilitate further discussion we introduce a $p\times 1$ column vector,
$|x\rangle$ of alternating 2’s and 0’s, for which the last element is a 2.

$|x\rangle=\left[\begin{array}[]{c}\vdots\\\ 0\\\ 2\\\ 0\\\
2\end{array}\right].$ (314)

From inspection of (314) and (36) we see that

$\hat{\mu}|x\rangle=\left[\begin{array}[]{c}\vdots\\\ 2\\\ 0\\\ 2\\\
0\end{array}\right].$ (315)

Therefore,

$(\hat{1}+\hat{\mu})|x\rangle=\left[\begin{array}[]{c}\vdots\\\ 2\\\ 2\\\ 2\\\
2\end{array}\right].$ (316)

The inverse $\hat{g}^{-1}$ of the statistical-weight matrix $\hat{g}$ of (164)
has the matrix elements

$(\kern-2.59995pt(l|\hat{g}^{-1}|l^{\prime})=g_{l}^{-1}\delta_{ll^{\prime}}.$
(317)

Using (317) with (300) and (316) we find

$2\hat{g}^{-1}|g\rangle=\left[\begin{array}[]{c}\vdots\\\ 2\\\ 2\\\ 2\\\
2\end{array}\right]=(\hat{1}+\hat{\mu})|x\rangle.$ (318)

Multiplying both sides of (318) by $(\hat{1}+\hat{\mu})^{-1}$, we find

$\displaystyle|x\rangle$ $\displaystyle=$ $\displaystyle
2(\hat{1}+\hat{\mu})^{-1}\hat{g}^{-1}|g\rangle$ (319) $\displaystyle=$
$\displaystyle 2[\hat{g}(\hat{1}+\hat{\mu})]^{-1}|g\rangle$ $\displaystyle=$
$\displaystyle 2\hat{h}^{-1}|g\rangle$ $\displaystyle=$
$\displaystyle|f\rangle.$

So the column vector $|f\rangle$ of (313) is the same as the column vector
$|x\rangle$ of (314), for which the elements alternate between 2 and 0, with
$x_{p-1}=2$.

In the eigenvalue equation of (312), the matrix $\hat{m}=|f\rangle\langle g|$
is the outer product between the column vector $|f\rangle$ and row vector
$\langle g|$. It is easy to find the eigenvectors $|\phi)$ and eigenvalues $m$
of outer products. There are $p-1$ degenerate eigenvalues, $m=0$,
corresponding to any set of $p-1$ independent vectors $|\phi)$ for which
$\langle g|\phi)=0$. These solutions are of no interest to us. The eigenvector
corresponding to the single, potentially non-zero eigenvalue can be chosen to
be

$|\phi\rangle=c|f\rangle.$ (320)

We will determine the constant $c$ below. Substituting (320) into the
eigenvalue equation of (312) we find

$c|f\rangle\langle g|f\rangle=c\,m|f\rangle.$ (321)

Eq. (321) implies that

$\langle g|f\rangle=m.$ (322)

Using (322) with the explicit form of $|f\rangle=|x\rangle$ of (314) and
$\langle g|$ of (300) we find that the value of the Lagrange multiplier $m$
must be

$\displaystyle m$ $\displaystyle=$ $\displaystyle
2(g_{p-1}+g_{p-3}+g_{p-5}+\cdots)$ (323) $\displaystyle=$ $\displaystyle
p(p+1).$

The sum on the first line of (323) only includes statistical weights $g_{l}$
for which $l\geq 0$. One can prove the last line of (323) by mathematical
induction. Letting $m_{p}$ be the value of $m$ for the phase function (298)
that is expanded on $2p$ Legendre polynomials, we see that (323) is true for
$p=1$ and $p=2$,

$\displaystyle m_{1}$ $\displaystyle=$ $\displaystyle 2(1)$ $\displaystyle
m_{2}$ $\displaystyle=$ $\displaystyle 2(3).$ (324)

If $m_{p}=p(p+1)$ for the integer $p$, then we can use (323) to write

$\displaystyle m_{p+2}$ $\displaystyle=$ $\displaystyle 2g_{p+1}+p(p+1)$ (325)
$\displaystyle=$ $\displaystyle 2(2p+3)+p^{2}+p$ $\displaystyle=$
$\displaystyle p^{2}+5p+6$ $\displaystyle=$ $\displaystyle(p+2)(p+3).$

This completes the proof that (323) is true for all positive integers $p$.

Substituting the coefficients $\phi_{l}=cf_{l}$ from (320) into (299) we find

$\displaystyle\phi(\mu)$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{p-1}\phi_{l}\,g_{l}\,P_{l}(\mu)$ (326)
$\displaystyle=$ $\displaystyle
2c\left[g_{p-1}P_{p-1}(\mu)+g_{p-3}P_{p-3}(\mu)+\cdots\right]$
$\displaystyle=$ $\displaystyle 2c\frac{dP_{p}}{d\mu}(\mu).$

An identity of Legendre polynomials [12] was used to go from the second to the
third line of (326). Substituting (326) into (298) we find

$\varpi^{\\{p\\}}(\mu)=4c^{2}(1+\mu)\left[\frac{dP_{p}}{d\mu}(\mu)\right]^{2}.$
(327)

Evaluating (327) for forward scattering ($\mu=1$) and using (311) with (323)
we find

$\varpi^{\\{p\\}}(1)=8c^{2}\left[\frac{dP_{p}}{d\mu}(1)\right]^{2}=p(p+1).$
(328)

The derivatives of the Legendre polynomials at the end point $\mu=1$ have the
values [12]

$\frac{dP_{p}}{d\mu}(1)=\frac{p(p+1)}{2}.$ (329)

Using (329) in (328) we find that the squared coefficient $c^{2}$ of (327)
must be

$c^{2}=\frac{1}{2p(p+1)}.$ (330)

Substituting (330) into (327) gives (141) and completes the proof that
$\varpi^{\\{p\\}}(\mu)$ is the phase function constructed from the first
$2p=2,4,6,\ldots$ Legendre polynomials that has the maximum possible forward-
scattering value.

## Acknowledgements

The Canadian Natural Science and Engineering Research Council provided
financial support of one of us.

## References

  * [1] C. Mobley, The Evolution of Radiative Transfer Theory.
https://www.oceanopticsbook.info/packages/iws_l2h/conversion/files/Mobley_EvolutionRTT_draft.pdf

  * [2] S. Chandrasekhar, Radiative Transfer, Dover, New York (1960).
  * [3] G. W. Petty A First Course in Atmospheric Radiation, Sundog Publishing, Madison, Wisconsin (2006).
  * [4] G. E. Thomas and K. Stamnes, Radiation Transfer in the Atmosphere and Ocean, Cambridge University Press, Cambridge, England (1999).
  * [5] A. D. Code, A Solution of a Time-Dependent Equation of Radiative Transfer, The Astrophysics Journal, 159, 1029 (1979).
  * [6] Reciprocal Lattice, http://www.phys.nthu.edu.tw/~spin/course/104F/solid%20state%20physics-Ch2-2-Tina.pdf
  * [7] W. A. van Wijngaarden and W. Happer, Dependence of Earth’s Thermal Radiation on Five Most Abundant Greenhouse Gases, http://arxiv.org/abs/2006.03098
  * [8] A. Schuster, Radiation Through a Foggy Atmosphere, The Astrophysical Journal, 21, p. 1 (1905).
  * [9] Gauss-Legendre quadrature,
https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature

  * [10] Magic angle, https://en.wikipedia.org/wiki/Magic_angle
  * [11] Kirchhoff’s laws of thermal radiation, https://en.wikipedia.org/wiki/Kirchhoff%27s_law_of_thermal_radiation
  * [12] Legendre Polynomials, https://en.wikipedia.org/wiki/Legendre_polynomials