\section{Introduction}

The success of LibriSpeech~\cite{povey2015} as a standard, freely available, Automatic Speech Recognition (ASR) benchmark is undeniable in the research community. LibriSpeech is English-only, and while benchmarks for other languages are available,  there are often low-scale or scattered around different places, and rarely available under an open license. In this paper, we revisit the work which has been done with LibriSpeech but in a multi-lingual manner and at a larger scale, introducing the Multilingual LibriSpeech (MLS) dataset. MLS includes 44.5K hours of English, and a total of 6K hours spread over 7 other languages. As for LibriSpeech, MLS is a read-speech dataset, which leverages LibriVox\footnote{\url{https://librivox.org}} audiobook data, most of which being based on the
Project Gutenberg\footnote{\url{http://www.gutenberg.org}} text data. LibriVox and Project Gutenberg data are released in the public domain, which allows us to release MLS freely to everyone.

In Section 3, we detail how we created the dataset, by (i)~training some acoustic models on in-house data, (ii)~generating pseudo-labels with these models, and (iii)~retrieving the original transcript by matching pseudo-labels to available book transcripts. Section~4 details the creation of train, development and test splits and limited supervision train sets. Section~5 introduces languages models we trained for each of language. These languages models are part of the MLS release. Section~6 covers some baseline ASR experiments.

\section{Related Work}

As for our work, LibriSpeech~\cite{povey2015} is derived from the LibriVox data, and is distributed under an open license. It ships with about $1000$ hours of labeled audio, obtained by leveraging alignments between textbooks and their read (audio) counterpart. In contrast to our work, it is only mono-lingual (English). A notable multi-lingual ASR dataset was built with the IARPA Babel Program~\cite{babel}. It collected data for 24 languages, mostly from conversational telephone speech. The dataset is however not released and under an open license, and focused on low-resource languages, with labeled data ranging between $25$ to $65$ hours per language. On the open license side, two important volunteer-supported multi-lingual speech gathering efforts are being conducted: (i)~VoxForge~\cite{Voxforge.org} which collected data for about 15 different languages, but remains low-scale (about $300$ hours in total). (ii)~CommonVoice~\cite{ardila2019common}, a more scalable solution, with more than 30 languages available, which keeps growing with $4500$ (validated) hours currently available. Other notatable multi-lingual datasets distributed under an open license are the M-AILABS~\cite{ailabs2019} and the CMU Wilderness~\cite{black2019} datasets. M-AILABS is a lower-scale version of our work, with 9 languages collected from LibriVox, for a total of about $1000$ hours available. The CMU Wilderness collects readings from the New Testament, with 700 different languages available.

\section{Data processing pipeline}
This section describes the major steps involved in preparing the MLS dataset.

\begin{table*}
\centering
 \caption{LibriVox audiobooks statistics for the top 15 languages (* - audiobooks with mix of multiple languages)} \label{tab:lv_data}
  \begin{tabular}{lrrr}
    \toprule
    \textbf{Language}      & \textbf{Hours}  & \textbf{Books}     & \textbf{Speakers}                \\
    \midrule
	English	& 71,506.79&12421&		4214 \\
	German	& 3,287.48& 593&		244 \\
	Dutch&		2,253.68& 206&	91 \\
	Spanish& 1,438.41&	285&		120 \\
	French&		1,333.35& 224&	114 \\
	Multilingual*& 516.82&	130	&	19 \\
	Portuguese& 	284.59	&	68&31 \\
	Italian	& 	279.43& 61&	28 \\
	Russian& 172.34&	 	44&	29 \\
	Latin& 138.93&		20&	16 \\
	Polish&		137.00& 25&	16 \\
	Church Slavonic&	136.42&		8&2 \\
	Hebrew	&125.72& 23	&	13 \\
	Japanese&97.67&	38&		24 \\
	Ancient Greek&69.77&	43&		8 \\
    \bottomrule
  \end{tabular}
\end{table*}

\subsection{Downloading audiobooks}
Table \ref{tab:lv_data} shows the LibriVox audiobooks data available for each language measured using LibriVox APIs~\footnote{\url{https://librivox.org/api/info}}. While English is the most dominant language, there is a significant amount of audio hours present in languages other than English, making this a valuable source for multilingual dataset preparation.

Based on the number of audiobook hours and the availability of the corresponding text sources of the audiobooks, we have selected English, German, Dutch, Spanish, French, Portuguese, Italian, Polish for the MLS dataset preparation. For downloading the LibriVox audiobooks in these languages, we have used data preparation tools available at Libri-Light~\footnote{\url{https://github.com/facebookresearch/libri-light}} open source library. All the audio data is downsampled from 48kHz to 16kHz for further processing.

\subsection{Audio segmentation}
Audio files downloaded from LibriVox have a huge variance in duration - from few minutes to hours.  Since acoustic model training is usually done on shorter utterances, we segment the audio files into 10-20 second segments. For this, we have used trained acoustic models in each of the languages using Time-Depth Separable (TDS) Convolutions\cite{Hannun2019SequencetoSequenceSR} with Auto-Segmentation Criterion~\cite{collobert2016wav2letter} (ASG) loss. We chose ASG criterion over Connectionist Temporal Classification~\cite{graves2006connectionist} (CTC) criterion since ASG doesn't exhibit delay in transcriptions compared to CTC~\cite{liptchinsky2017letterbased}. For each language, we train models on in-house datasets consisting of videos publicly shared by users. We use only audio part of the videos and the data is completely de-identified.

Our segmentation process is shown in Figure \ref{fig:segment} and consists of two steps. First, we run inference on the audio and generate viterbi token sequence along with their timestamps. Since the audio files can be very long, we have used wav2letter@anywhere framework~\cite{pratap2020scaling} to perform the inference in a streaming fashion. Second, we select the longest silence duration within 10 sec to 20 sec range from the start of an audio and split the audio at the mid point of the silence chunk to create a segment. If no silence frames are found between 10 sec to 20 sec from the starting point, we split the audio at 20 sec mark. Once we generate a segment, we consider the end point of previous segment as the current starting point and repeat the process again till we reach the end of audio.  This process guarantees that all the segments are between 10sec and 20sec. A minimum segment duration of 10 sec is kept so that the segments have sufficient number of words spoken which helps with better transcript retrieval (described in Section ~\ref{trans_ret}).

\begin{figure}
\centering
   \includegraphics[width=\linewidth]{images/segment.pdf}
  \caption{Audio segmentation procedure to split long audio files into shorter segments}
  \label{fig:segment}
\end{figure}

\subsection{Pseudo label generation}
\label{sec:pl}
We generate pseudo labels for the segmented audio samples by performing a beam-search decoding with a 4-gram language model on the same acoustic models used for audio segmentation. The language models are trained on the data used for training the acoustic models. For English, however, we use pre-trained model from~\cite{synnaeve2019endtoend} which uses TDS encoder and CTC loss on LibriSpeech~\cite{panayotov2015librispeech} and pseudo labels from LibriVox .

\subsection{Downloading text sources for audiobook data}
To generate the labels for the audio segments derived from audiobooks, we would need to have the original textbook from which the speaker read the audiobook. For English, we found that $\approx$60K hours of audiobooks is read from four major website domains - \url{gutenberg.org}, \url{archive.org}, \url{ccel.org} and \url{hathitrust.org}. We have written parsers to automatically extract the text for each of these domains and downloaded the text sources for all the audiobooks in English.

For other languages, however, we found it more challenging to automatically extract text for the audiobooks because: 1. the diversity of domains is large making it difficult to write parsers for each and every domain 2. some of the links were invalid or redirected to an incorrect page. So, we have also incorporated some manual approaches in our process to cover the audiobook text sources as much as possible. Depending on the language and text source domain, we copied the text data directly from the  browser or extracted text from .pdf/.epub books using pdftotext\footnote{\url{https://pypi.org/project/pdftotext/}} tool, or wrote HTML parsers to retrieve text data for popular domains in a language.

For the audiobooks with invalid links for text sources, we have also manually searched online to find alternate sources where the book is available. For example, all the text sources from the \url{spiegel.de} domain, which accounts of 1/3rd of German audio data, were being redirected to an invalid page. However, we were able to find the alternate text sources from websites like \url{projekt-gutenberg.de}, \url{zeno.org} for most of these unavailable books from \url{spiegel.de}.

\subsubsection{Text Normalization}
\label{tnorm}
Given that the text data is taken from many different online sources, we perform careful text normalization to make the data easy to work with for later stages of the data processing pipeline. For normalizing the text, we first perform NFKC normalization and remove all the unwanted characters like punctuations, subscript/superscripts, emojis, escape symbols etc.  We also remove the hyphens used for ``end-of-line hyphenation"~\footnote{\url{https://practicaltypography.com/hyphenation.html}} and join the parts of words into a single word. We then prepare a list of valid unicode characters based on the language’s orthography and filter characters outside this range.

\subsection{Transcript retrieval}
\label{trans_ret}

The transcript retrieval process involves finding the true target label for the audio segments from the source text of audio. Our procedure closely follows the method described in~\cite{manohar2017jhu} with few modifications.

We first split the source text into multiple overlapping documents of 1250 words each and striding by 1000 words. We retrieve the documents which best matches with the pseudo label for the audio segments using term-frequency inverse document-frequency (TF-IDF) similarity score on bigrams. We then perform a Smith-Waterman alignment~\cite{SMITH1981195} to find the best matching sub-sequence of words. We have used a matching score of 2 and substitution, insertion, deletion score of -1 for the alignment algorithm.

After the above alignment procedure, we generate a candidate target label for each audio segment, which corresponds to the best match of the pseudo label in the source text of audiobook. We filter out all the candidate transcripts generated from the matching algorithm above, if the WER between the candidate transcript and pseudo label generated is \textgreater 40\%.

\subsubsection{Post processing of numbers, hyphens, apostrophe }

\begin{figure}
\begin{minipage}{.48\textwidth}
  \centering
    \caption{Example implementation of number replacement procedure. Numbers from book text are replaced with the aligned word from pseudo label.}
     \includegraphics[width=0.95\linewidth]{images/numbers.pdf}
  \label{fig:numbers}
\end{minipage}%
\hfill
\begin{minipage}{.48\textwidth}
  \centering
    \caption{Heuristics used for dealing with hyphens and apostrophes. Words marked in bold are replaced with the words present below them.}
     \includegraphics[width=0.95\linewidth]{images/hyphen.pdf}
  \label{fig:hyphens}
\end{minipage}
\end{figure}

Numbers in the text present a unique challenge in the preparation of the dataset. While a library like num2words\footnote{\url{https://pypi.org/project/num2words/}} can be used to convert numbers to words and supports multiple languages, we found that using such tools won't solve the problem in a generic way. For example, a number ``401" could be spelled as ``four hundred and one"  or  ``four-o-one" depending on the context and it is very hard to come up with automated rules to pick one over the other. Further, some of the numbers present in the text could be page numbers and are not read by the speaker. To solve the issue, we align the matched text from the book with the pseudo label and replace the numbers in the matched text with the corresponding aligned words from the pseudo label. While this solution is not perfect as the pseudo label may be not be always have the correct transcription of the audio, we found this to be a reasonably good solution to solve the issue. Figure \ref{fig:numbers} shows an example implementation of this procedure.

On the other hand, dealing with hyphens and apostrophes is challenging because we found these character to appear in the text in many places where we do not expect them to be for an audio transcript. Figure \ref{fig:hyphens} shows some examples of these issues from the books we processed and the heuristics we used to replace such word forms. To determine is a word is rare, we count of number of distinct books the word appears and check if it is less than a chosen threshold.

\section{Putting it together}
\subsection{Creating train, development and test splits}

\begin{table*}[]
\centering
\caption{Statistics of Train/Dev/Test partitions of each language. Below lists for each partition: the total duration in hours (left), number of speakers in each gender (middle) and duration per gender in dev and test sets (right). \label{tab:speaker_stat}}
\begin{tabular}{@{}cccccccccccccc@{}}
\toprule
\multirow{3}{*}{\textbf{Language}} & \multicolumn{3}{c}{\textbf{Duration (hrs)}}                                          & \multicolumn{6}{c}{\textbf{\# Speakers} }                                                & \multicolumn{4}{c}{\textbf{\# Hours / Gender}}             \\ \cmidrule(lr){2-4} \cmidrule(lr){5-10} \cmidrule(lr){11-14}
                          & \multirow{2}{*}{train} & \multirow{2}{*}{dev} & \multirow{2}{*}{test} & \multicolumn{2}{c}{train} & \multicolumn{2}{c}{dev} & \multicolumn{2}{c}{test} & \multicolumn{2}{c}{dev} & \multicolumn{2}{c}{test} \\ \cmidrule(lr){5-6} \cmidrule(lr){7-8} \cmidrule(lr){9-10} \cmidrule(lr){11-12} \cmidrule(lr){13-14}
                          &                        &                      &                       & M           & F           & M          & F          & M           & F          & M          & F          & M           & F          \\
\midrule
English    & 44,659.74 & 15.75 & 15.55 & 2742 & 2748 & 21 & 21 & 21 & 21 & 7.76 & 7.99 & 7.62 & 7.93 \\
German     & 1,966.51  & 14.28 & 14.29 & 81   & 95   & 15 & 15 & 15 & 15 & 7.06 & 7.22 & 7.00    & 7.29 \\
Dutch      & 1,554.24  & 12.76 & 12.76 & 9    & 31   & 3  & 3  & 3  & 3  & 6.44 & 6.32 & 6.72 & 6.04 \\
French     & 1,076.58  & 10.07 & 10.07 & 62   & 80   & 9  & 9  & 9  & 9  & 5.13 & 4.94 & 5.04 & 5.02 \\
Spanish    & 917.68    & 9.99  & 10    & 36   & 50   & 10 & 10 & 10 & 10 & 4.91 & 5.08 & 4.78 & 5.23 \\
Italian    & 247.38    & 5.18  & 5.27  & 22   & 43   & 5  & 5  & 5  & 5  & 2.5  & 2.68 & 2.38 & 2.90  \\
Portuguese & 160.96    & 3.64  & 3.74  & 26   & 16   & 5  & 5  & 5  & 5  & 1.84 & 1.81 & 1.83 & 1.90  \\
Polish     & 103.65    & 2.08  & 2.14  & 6    & 5    & 2  & 2  & 2  & 2  & 1.12 & 0.95 & 1.09 & 1.05 \\
\bottomrule
\end{tabular}
\end{table*}

We have used the following principles when splitting the dataset into train, development and test sets - 1) there is no speaker overlap between the training, development and test sets, 2) speakers are balanced in gender and duration in development and test sets and 3) there are sufficient audio hours and speakers assigned into development and test sets to be able to validate ASR model performance.

First, we select the list of all the books available for a language. We remove all the books with corrupted meta data, such as missing title or information about speakers and authors. Then, to ensure that each recording is unambiguously attributable to a single speaker, we also remove audios with multiple speakers, for example, “Dramatic Reading”, which include predominantly multi-reader audio chapters. In addition, we only keep the latest version for books sharing the same authors and title, but different versions. Only the speakers reading those valid books are considered in the transcript retrieval process from section \ref{trans_ret}.

Second, we label the gender of all the speakers with a gender classifier. The classifier is a SVM \cite{svm} with RBF kernel trained on 1172 speakers from \textit{train-clean-100} and \textit{train-clean-360} subsets of LibriSpeech . In particular, it consumes 40-dimensional log-filterbank features averaged over time as input features. The test accuracy on the 146 speakers from the joint development and test sets of LibriSpeech is 95\%. We use the same gender classifier for other languages as well. We manually checked the quality of this classifier on Dutch and Polish and its accuracy is 96\% and 94\% respectively.

Finally, the dataset is split into training, development and test sets as following. We computed the total duration each speaker spends in reading the valid books, and order the speakers by this duration. Speakers with duration shorter than a threshold are assigned into the training set. Then, from the rest of the speakers, we select a series of speakers with the shortest duration in each gender equally. The selected speakers are then equally split into development and test sets. All the remaining speakers are assigned to training set again. Finally, after the speakers are assigned to their given partitions, we simply attach the samples recorded by them into different partitions accordingly. To avoid high speaker imbalance in development and test sets, we further truncate speakers with high duration by sampling their recordings up to an upper-bound. We also make sure that each chapter of each valid book only appears in one partition. The detailed statistics can be found in Table \ref{tab:speaker_stat}.

For English, aside from the same procedure described above, we additionally make it exclusive to the previous LibriSpeech dataset. As shown in Figure \ref{fig:partitions}, there is speaker and book overlap between the two train sets, but the development and test sets are completely independent from one another. This makes the 2 training sets interchangeable and complementary. Note that there is book overlap across partitions, because following LS preparation, different chapters of a given book may be read by different speakers. In addition, the MLS development and test sets are chosen to be harder than LibriSpeech dev-/test-other paritions. Specifically, by comparing the pseudo labels generated as section \ref{sec:pl} with the final transcriptions, speakers with WER higher than 80\% of the ones in dev-other sets are picked into MLS dev/test sets.

\begin{figure}
\centering
\begin{minipage}{.48\textwidth}
  \centering
    \caption{Speaker (fully filled), chapter (gridline filled) and book (empty) overlap of MLS and LibriSpeech (LS) partitions for English. The figure is not to scale.}
   \includegraphics[width=0.95\linewidth]{images/partitions.pdf}
  \label{fig:partitions}
\end{minipage}%
\hfill
\begin{minipage}{.48\textwidth}
  \centering
    \caption{Violin plots of  audio segments duration in the training data for different languages}
    \includegraphics[trim={1.1cm 0.5cm 1.8cm 0.5cm},clip, width=\linewidth]{images/violin_plot.png}

  \label{fig:violin}
\end{minipage}
\end{figure}

\subsection{Limited supervision train data}
We also provide limited supervision dataset for each language to create a standard benchmark for low resource training with MLS. Similar to \cite{9052942}, we provide a 10 hour set, a 1 hour set, and six 10-minute sets (the six 10-minute sets together make up the 1h set, and the 1h set is included in the 10h set). For preparing the six 10-minute sets, we first sample upto 15 male and 15 female speakers from train set for each language. From this set of speakers, we select 3 male speakers and 3 female speakers at random and sample 5 minutes of audio for each gender making sure we don't select a sample from a previous 10 min train set. For preparing the remaining 9 hour set, we consider all the speakers selected initially and sample 4.5 hours of audio per gender that is not present in the six 10-minute sets created before.

\subsection{High quality development and test sets}
Using a high quality development and test sets are critical for measuring ASR system performance and making decisions based on the test results. In order to make sure the quality of MLS development and test sets is high, we perform human evaluation of the transcripts in these subsets. The human transcribers listen to the audio file while looking the transcript and are asked to correct the transcripts if they find any mistakes. To ensure that we do not add any potential bias into the dataset from human transcribing, we ask the transcribers to only correct ``clear and obvious" errors.

Figure \ref{fig:human_corrections} shows some of the examples of corrections using this workflow and Table \ref{tab:human_wer} shows the WER between generated transcript and human rated transcript.

\begin{table}
\centering
\begin{minipage}{.48\textwidth}
  \centering
    \captionof{figure}{Example corrections to generated transcriptions made as part of human verification workflow. GT - Generated Transcript, HRT - Human Rated Transcript}
      \begin{tcolorbox}[colback=gray!5!white,colframe=white!20!black,boxrule=0.5pt,arc=0pt, boxsep=2pt,left=2pt,right=2pt,top=2pt,bottom=2pt]
   \textit{\# Diactric correction} \\
   GT:  menesteroso para matar á los de recto proceder \\
   HRT: menesteroso para matar a los de recto proceder \\

   \textit{\# Missing words} \\
   GT: cabeça escandecida amor que vive\\
   HRT: cabeça escandecida amor que vive e brilhe \\

   \textit{\# Missing Apostrophe} \\
   GT: teseo combatter co doppi petti e degli ebrei\\
   HRT: teseo combatter co' doppi petti e degli ebrei
\end{tcolorbox}

  \label{fig:human_corrections}
\end{minipage}%
\hfill
\begin{minipage}{.48\textwidth}
\centering
\caption{WER between generated transcript (GT) and human rated transcript (HRT). GT is considered as hypothesis and HRT as reference for measuring WER.}
\begin{tabular}{ccc}
    \toprule
        \multirow{2}{*}{\textbf{Language}} & \multicolumn{2}{c}{\textbf{Error Rate}} \\
 \cmidrule(lr){2-3}
        & Dev & Test \\
    \midrule
     English&	4.55 & 	4.54 \\
German & 	3.6 & 	3.9 \\
Dutch & 	7.24 & 	8.61 \\
French & 	3.33 & 	4.45 \\
Spanish & 	1.93 & 	3.07 \\
Italian & 	6.8 & 	5.37 \\
Portuguese & 	12.64 & 	12.07 \\
Polish & 	5.67 & 	5.12 \\
    \bottomrule
\end{tabular}
\label{tab:human_wer}
\end{minipage}
\end{table}

\subsection{Dataset statistics}
Table \ref{tab:speaker_stat} shows the amount of train, valid and test data  along with the gender distribution of speakers for the 8 languages that we processed.  The training data is lower than the data we presented in Table \ref{tab:lv_data} because of 1) we may not have downloaded the text source for the audiobook, 2) the candidate label is filtered during the transcript retrieval stage or 3) filtered while creating valid/test splits.

 Figure~\ref{fig:violin} shows the violin plots of the duration of audio segments in training data for each language. We can see that all the segments are within 10sec to 20 sec range and all the sizes are almost evenly distributed inside the range.

\section{Language models}
\label{sec:lm}
We also release language model training data and pre-built language models along with the MLS dataset. The source text for these language models is taken from Project Gutenberg\footnote{http://gutenberg.org/} books and the text is carefully normalized similar to the steps mentioned in Section \ref{tnorm}.

We have also carefully filtered certain books to avoid any overlap with the development and test sets, following the procedure described in \cite{panayotov2015librispeech}. Specifically, any book whose title has edit distance (over words) less than $2$ with any of the titles of books in development and test sets are removed in their entirety. Then, an inverted index of all 5-grams, with stop words removed, is built for the books in the test and development sets. All candidate books are then checked against this index and any book with more than one percent of the 5-grams which appear in it is removed as well.

We have trained 3-gram and 5-gram language models (LM) for all the languages in our dataset using the KenLM toolkit~\cite{heafield2011kenlm}. We release all the language model training data and pre-built language models with the dataset. The detailed information of language model corpus for each language can be found in Table \ref{tab:lm_corpus}. The  out-of-vocabulary (OOV) rate and their perplexities (excluding OOV words) on the transcriptions in development sets are listed in Table \ref{tab:lm_ppl}.

\begin{table*}[]
\centering
\caption{Language model corpus statistics. Number of books before and after filtering (col 2 - 3), number of words (col 4 - 5) and number of sentences (col 6). \label{tab:lm_corpus}}
\begin{tabular}{@{}cccrrr@{}}
\toprule
\textbf{Language} & \textbf{\# Books} & \textbf{\# Books (after filtering)}  & \textbf{\# Words} & \textbf{\# Words (Unique)} & \textbf{\# Sentences} \\
\midrule
English           & 47677            & 36866 (77\%)                 & 2.38B         & 4.12M                     & 133.26M             \\
German            & 1757             & 1323 (75\%)                 & 67.06M           & 1.24M                     & 3.81M               \\
Dutch             & 796              & 629 (79\%)                 & 36.26M           & 559.6K                      & 1.81M               \\
French            & 3008             & 2424 (81\%)                 & 164.1M          & 911.22K                      & 8.34M               \\
Spanish           & 622              & 479 (77\%)                 & 32.41M           & 489.94K                      & 1.61M               \\
Italian           & 814              & 640 (79\%)                 & 42.28M           & 700.63K                      & 2.14M               \\
Portuguese        & 554              & 460 (83\%)                 & 13.84M           & 383.2K                      & 1.24M               \\
Polish            & 31               & 22 (71\%)                 & 492.32K             & 67.1K                       & 28.86K             \\
\bottomrule
\end{tabular}
\end{table*}

\begin{table*}[]
\centering
\caption{The performance of the language models of each language.  Out-of-vocabulary (OOV) rates on dev sets (col 2) and perplexity on dev sets for 3-gram and 5-gram LMs (col 3). \label{tab:lm_ppl}}
\begin{tabular}{@{}cccc@{}}
\toprule
 & & \multicolumn{2}{c}{\textbf{LM Perplexity}} \\
\cmidrule(lr){3-4}
\multirow{-2}{*}{\textbf{Language}} & \multirow{-2}{*}{\textbf{OOV rate}} & 3-gram  & 5-gram \\
\midrule
English     & 0.18\% & 232.11 &   190.76            \\
German      & 1.22\% & 592.22  &   570.84            \\
Dutch       & 1.50\% & 461.22 &   439.62            \\
French      & 0.42\% & 261.65 &   238.47            \\
Spanish     & 0.92\% & 376.81       &    358.09     \\
Italian     & 1.72\% &  733.10      &     632.53    \\
Portuguese  & 1.81\% & 983.71       &    949.96   \\
Polish      & 13.39\%  &  2450.23      &   2442.12  \\
\bottomrule
\end{tabular}
\end{table*}

\section{Experiments and results}

\subsection{Monolingual Baselines}

All the baseline models are trained using wav2letter++ \cite{pratap2018} framework. The  network architecture of acoustic model closely follows \cite{likhomanenko2020rethinking}: the encoder of our AMs is composed of a convolutional frontend (1-D convolution with kernel-width 7 and stride 3 followed by GLU activation) followed by 36 4-heads Transformer blocks\cite{vaswani2017attention}. The self-attention dimension is 768 and the feed-forward network (FFN) dimension is 3072 in each Transformer block. The output of the encoder is followed by a linear layer to the output classes. We use dropout after the convolution layer. For all Transformer layers, we use dropout on the self-attention and on the FFN, and layer drop~\cite{fan2020}, dropping entire layers at the FFN level. Dropout and layer dropout values are tuned for each model separately.

We use the Adagrad optimizer~\cite{duchi2011adaptive} and decay learning rate by a factor of 2 each time the WER reaches a plateau on the validation sets.
The AMs take 80-channel log-mel filterbanks as input and are trained end-to-end with Connectionist Temporal Classification (CTC) loss~\cite{graves2006connectionist}. SpecAugment\cite{Park_2019} is used for data augmentation in training: there are two frequency masks, and ten time masks with maximum time mask ratio of p = 0.1; frequency and time mask parameters are tuned separately for each model; time warping is not used. The maximum frequency bands masked by one frequency mask is 30, and the maximum frames masked by the time mask is 30, too. We use the set of graphemes from the corresponding training set of the language as the token set. The viterbi WER on development and test for all the languages are shown in Table \ref{tab:baslines}.

To further improve the WER, we use beam-search decoding in wav2letter++ without using a LM (ZeroLM) and using 5-gram language model mentioned in Section \ref{sec:lm}. The decoder hyper-parameters are tuned on the development set and the results are shown in Table \ref{tab:baslines}. We can see that WER improves when decoding with a 5-gram LM for all the languages except Polish. For Polish, the OOV rate is 13\% as shown in Table \ref{tab:lm_ppl} which negatively impacts the WER as we are performing LM decoding with constrained lexicon. All the pretrained models and the recipes to train the acoustic model and decoding them will be made available at \url{https://github.com/facebookresearch/wav2letter}.

\begin{table*}[]
\centering
\caption{Monolingual Baseline WER with different decoding strategies. \label{tab:baslines}}

\begin{tabular}{@{}ccccccc@{}}
\toprule
                                    & \multicolumn{2}{c}{\textbf{Viterbi}} & \multicolumn{2}{c}{\textbf{Zero LM}} & \multicolumn{2}{c}{\textbf{5-gram LM}} \\ \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7}

\multirow{-2}{*}{\textbf{Language}} & dev               & test             & dev                          & test  & dev                 & test                \\
\midrule
English                             & 6.01              & 6.99             & 5.89                         & 6.76  & 5.05                & 5.88                \\
German                              & 5.55              & 6.93             & 5.80  & 7.10   & 5.45                & 6.49                \\
Dutch                               & 17.00                & 13.18            & 16.45                        & 13.09 & 14.42               & 12.02               \\
French                              & 7.79              & 6.88             & 7.43                         & 6.58  & 6.58                & 5.58                \\
Spanish                             & 5.94              & 6.90              & 5.90                          & 6.68  & 5.27                & 6.07                \\
Italian                             & 14.55             & 12.35            & 14.01                        & 11.78 & 11.96               & 10.54               \\
Portuguese                          & 18.62             & 21.70             & 17.22                        & 20.52 & 16.35               & 19.49               \\
Polish                              & 19.25             & 19.40             & 18.73                        & 21.66 & 17.64               & 20.39              \\
\bottomrule
\end{tabular}
\end{table*}

\subsection{Comparison with LibriSpeech}

\begin{table*}[]
\centering
\caption{Comparing model performance on Librispeech dev/test sets using LibriSpeech and MLS (different amount of training data, same model setup). \label{tab:ls_wer}}
\begin{tabular}{@{}cccccccc@{}}
\toprule
\multirow{2}{*}{\textbf{Training Dataset}} &
\multirow{2}{*}{\textbf{Training Hours}}&
\multicolumn{2}{c}{\textbf{Viterbi}} & \multicolumn{2}{c}{\textbf{4-gram LM}} & \multicolumn{2}{c}{\textbf{Transformer LM}} \\ \cmidrule(lr){3-4} \cmidrule(lr){5-6} \cmidrule(lr){7-8}
                                      &  & clean             & other            & clean                    & other                   & clean                     & other                     \\
\midrule
LibriSpech  \cite{likhomanenko2020rethinking}      &     960                   & 2.8 / 2.8           & 7.1 / 7.1          & 2.0 / 2.5                  & 5.3 / 5.6                 & 1.5 / 2.1                   & 4.3 / 4.7                   \\
MLS                   &  44.5K         & 2.94 / 3.09       & 5.39 / 5.50       & 2.06 /2.33                & 4.01 / 4.38               & 1.83 / 2.11                 & 3.54 / 3.97                \\
\bottomrule
\end{tabular}
\end{table*}

Table \ref{tab:ls_wer} compares the WER performance on LibriSpeech dev and test sets using the same model architecture. We use the 4-gram LM provided with LibriSpeech for decoding both models and the transformer LM used for rescoring is trained on the LM corpus and training data provided with LibriSpeech. It can be seen that using MLS dataset for English and thereby using large amount of supervised training data improves the performance on LibriSpeech dev and test sets.

\section{Conclusions}

We have presented the Multilingual LibriSpeech dataset, a large scale multilingual speech dataset with 36.5K hours of training data spread over 8 languages. We believe this dataset will promote open research in large-scale training of ASR systems and in multilingual ASR. This dataset can also be used for Text-to-Speech (TTS) research by extending the LibriTTS~\cite{48008} dataset, and by creating a larger and multilingual version for TTS Research.
\section{Acknowledgements}
We would like to thank Steven Garan for help in data preparation and text normalization and Mark Chou for helping with setting up the workflow for transcription verification.

\newpage

\bibliographystyle{IEEEtran}


\section{Introduction}
\label{sec:1}

Given a locally finite set $A \subseteq \Rspace^d$ and a line segment,
the \emph{Voronoi path} of the line segment is the dual of the Voronoi tessellation
of $A$ intersected with the segment.
In other words, it consists of all Delaunay edges dual to Voronoi cells
of dimension $d-1$ crossed by the line segment.
We generalize it to the \emph{Voronoi scape} of $A$ and a $p$-dimensional set
$\Omega \subseteq \Rspace^d$, which contains a $p$-cell in the Delaunay mosaic
of $A$ with multiplicity $\mu$, if $\Omega$ intersects the corresponding
Voronoi $(d-p)$-cell in $\mu$ points.
We are interested in the \emph{distortion}, which is the ratio of
the $p$-dimensional volume of the Voronoi scape over the volume of $\Omega$.

Considering the Voronoi tessellation of a stationary Poisson point process
and a line segment in $\Rspace^2$,
\cite{BTZ00} proves that the expected distortion is $\frac{4}{\pi}$.
Extending this work to $d > 2$ dimensions,
\cite{dCDe18} proves that the expected distortion is $\sqrt{2d/\pi} + O(1/\sqrt{d})$.
We remove the ambiguity in this answer by proving that
the expected distortion in $\Rspace^d$ is $d!!/(d-1)!!$, if $d$ is odd,
and $2 d!! / (\pi (d-1)!!)$, if $d$ is even, in which $!!$ is the double factorial.
Furthermore, we generalize the result to any dimension
and prove that for $p$-dimensional sets the expected distortion
is the binomial coefficient $\tbinom{d/2}{p/2}$, in which non-integer parameters
are understood in the way the Gamma function extends the factorial:
\begin{align}
    \DConstant{p}{d}  &=  \binom{d/2}{p/2}
      =  \frac{\Gamma(\frac{d}{2} + 1)}{ \Gamma(\frac{p}{2} + 1)
         \, \Gamma(\frac{d-p}{2} + 1)} \renewcommand{\arraystretch}{1.5}
      =  \left\{ \begin{array}{ll}
           \frac{d!!}{p!! \, (d-p)!!} \frac{2}{\pi}
             &  \mbox{\rm if $d$ is even and $p$ is odd} , \\
           \frac{d!!}{p!! \, (d-p)!!}  &  \mbox{\rm otherwise}.
       \end{array} \right.
    \label{eqn:DConstant}
\end{align}
The binomial interpretation also provides the asymptotics for $\DConstant{p}{d}$.
Table \ref{tbl:distortion} shows the values for small dimensions.
More precisely, we prove that \eqref{eqn:DConstant} is
the \emph{average distortion} for sufficiently regular $p$-dimensional sets
and Voronoi tessellations.
The claim for stationary Poisson point processes follows because
they are invariant under rotations and translations.
The proof is based on a decomposition of $\Rspace^d \times \Grass{p}{d}$
that is related to the \emph{mixed complex} \cite{Ede99}.
As a byproduct, we get an expression for the volumes of the cells
in the mixed complex; see Corollary \ref{cor:mixed_volumes}.
{\renewcommand{\arraystretch}{1.6}
\begin{table}[!hbt
  \centering \footnotesize
  \begin{tabular}{c||cccccccccc}
  & $d = 1$ & $2$ & $3$ & $4$ & $5$ & $6$ & $7$ & $8$ & $9$ & $10$ \\
 \hline \hline
 $p = 1$ & $1$ & $\frac{4}{\pi }$ & $\frac{3}{2}$ & $\frac{16}{3 \pi }$ & $\frac{15}{8}$ & $\frac{32}{5 \pi }$ & $\frac{35}{16}$ & $\frac{256}{35 \pi }$ & $\frac{315}{128}$ & $\frac{512}{63 \pi }$ \\
 $2$ &  & $1$ & $\frac{3}{2}$ & $2$ & $\frac{5}{2}$ & $3$ & $\frac{7}{2}$ & $4$ & $\frac{9}{2}$ & $5$ \\
 $3$ &  &  & $1$ & $\frac{16}{3 \pi }$ & $\frac{5}{2}$ & $\frac{32}{3 \pi }$ & $\frac{35}{8}$ & $\frac{256}{15 \pi }$ & $\frac{105}{16}$ & $\frac{512}{21 \pi }$ \\
 $4$ &  &  &  & $1$ & $\frac{15}{8}$ & $3$ & $\frac{35}{8}$ & $6$ & $\frac{63}{8}$ & $10$ \\
 $5$ &  &  &  &  & $1$ & $\frac{32}{5 \pi }$ & $\frac{7}{2}$ & $\frac{256}{15 \pi }$ & $\frac{63}{8}$ & $\frac{512}{15 \pi }$ \\
 $6$ &  &  &  &  &  & $1$ & $\frac{35}{16}$ & $4$ & $\frac{105}{16}$ & $10$ \\
 $7$ &  &  &  &  &  &  & $1$ & $\frac{256}{35 \pi }$ & $\frac{9}{2}$ & $\frac{512}{21 \pi }$ \\
 $8$ &  &  &  &  &  &  &  & $1$ & $\frac{315}{128}$ & $5$ \\
 $9$ &  &  &  &  &  &  &  &  & $1$ & $\frac{512}{63 \pi }$ \\
 $10$ &  &  &  &  &  &  &  &  &  & $1$
\end{tabular}
  \caption{The average resp.\ expected distortion in small dimensions.
           Note that even rows and columns form the Pascal triangle.}
  \label{tbl:distortion}
\end{table}}

\medskip \noindent \textbf{Outline.}
Section~\ref{sec:2} prepares the proof of our main result by computing
the first and second moments of the $p$-dimensional volume of the projection
of a unit $p$-cube in $\Rspace^d$.
Section~\ref{sec:3} studies the space of point-direction pairs.
Section~\ref{sec:4} introduces a mild regularity condition for Voronoi tessellations.
Section~\ref{sec:5} computes the volume of the cells in the mixed complex.
Section~\ref{sec:6} proves that $\DConstant{p}{d}$ is the average distortion
for $p$-dimensional shapes in $\Rspace^d$, and the expected distortion
if the tessellation is of a stationary Poisson point process.
Section~\ref{sec:8} concludes the paper.

\section{Random Projections}
\label{sec:2}

We need some preliminary computations.
Let $\Grass{p}{d}$ be the \emph{linear Grassmannian manifold},
whose points are the $p$-planes that pass through the origin of $\Rspace^d$.
Given a $p$-dimensional unit cube, $E \subseteq \Rspace^d$,
and a $p$-plane, $L \in \Grass{p}{d}$,
we write $\proj{E}{L}$ for the projection of the cube onto the plane,
and $\Vol{p}{\proj{E}{L}}{}$ for its $p$-dimensional volume.
The \emph{$j$-th projection moment} is the average $j$-th power
of the volume of the projection.
We express this moment as an integral over the Grassmannian
equipped with the uniform probability measure in \eqref{eqn:momentDef1}
and convert it to two equivalent expressions involving the angle
to a fixed plane in \eqref{eqn:momentDef2} and \eqref{eqn:momentDef3}:
\begin{align}
	\moment{p}{d}{j}  &=  \int\nolimits_{L \in \Grass{p}{d}} \VolPower{p}{\proj{E}{L}}{j} \diff L ,
	\label{eqn:momentDef1} \\
	&=  \int\nolimits_{L \in \Grass{p}{d}} \cos^j \Angle{L}{L_0} \diff L
	\label{eqn:momentDef2} \\
	&=  \int\nolimits_{F \in \Stiefel{p}{d}} \VolPower{p}{\proj{F}{L_0}}{j} \diff F .
	\label{eqn:momentDef3}
\end{align}
To explain \eqref{eqn:momentDef2} and \eqref{eqn:momentDef3},
we fix the plane $L_0 \in \Grass{p}{d}$ containing $E$.
The \emph{angle} between two $p$-planes, $\Angle{L}{L_0} \in [0, \frac{\pi}{2} ]$,
is defined as the arc-cosine of the ratio of $\Vol{p}{\proj{B}{L}}{}$
over $\Vol{p}{B}{}$ for any compact set with non-empty interior, $B \subseteq L_0$.
The angle is symmetric, so we can instead consider the integrand
in \eqref{eqn:momentDef2} as the projection of a unit $p$-cube
in a random $p$-plane onto $L_0$.
Formally, we write $\Stiefel{p}{d}$ for the \emph{Stiefel manifold}
of orthonormal $p$-frames in $\Rspace^d$,
we identify a frame with the unit $p$-cube it spans,
and we integrate using the uniform probability measure of $\Stiefel{p}{d}$
to arrive at \eqref{eqn:momentDef3}.

\medskip
By construction, the $0$-th projection moment is equal to $1$,
independent of $p$ and $d$.
We compute the $1$-st and $2$-nd projection moments, which curiously
both have intuitive geometric interpretations.
\begin{lemma}[Projection Moments]
  \label{lem:projection_moments}
  Let $d \geq 0$ and $0 \leq p \leq d$.
  Then
  \begin{align}
    \moment{p}{d}{1}  &=  \frac{\Gamma (\frac{p+1}{2}) \; \Gamma (\frac{d-p+1}{2})}
	{\Gamma (\frac{1}{2}) \; \Gamma (\frac{d+1}{2})} ,
    \label{eqn:moment1} \\
    \moment{p}{d}{2}  &=  1 / \binom{d}{p} =  \frac{p! \, (d-p)!}{d!}.
    \label{eqn:moment2}
  \end{align}
\end{lemma}
\begin{proof}
  The $1$-st projection moment appears in the classic Crofton formula
  of integral geometry, which says that the volume of a convex body
  is proportional to the average volume of its projections.
  The constant of proportionality given in \eqref{eqn:moment1} can be found
  in \cite[Formula (5.8)]{ScWe08}.

  We use \eqref{eqn:momentDef3} together with a generalization of
  the Pythagorean theorem to compute the $2$-nd moment.
  By Pythagoras, the squared length of a line segment is the sum
  of squared lengths of its projections onto the coordinate axes.
  The Cauchy--Binet formula \cite[\S4.6]{BrWi89} can be used to generalize
  this to the squared volume of a $p$-dimensional parallelepiped in $\Rspace^d$.
  Let $P$ be such a parallelepiped, and write $P_i$ for its projection
  onto the $i$-th coordinate $p$-plane (in which the numbering is arbitrary).
  There are $\binom{d}{p}$ coordinate $p$-planes, and the Cauchy--Binet formula asserts
  \begin{align}
    \VolPower{p}{P}{2}  &=  \sum\nolimits_{i=1}^{\binom{d}{p}} \VolPower{p}{P_i}{2} .
    \label{eqn:CauchyBinet}
  \end{align}
  Letting $P = F \in \Stiefel{p}{d}$ be the uniformly random unit $p$-cube,
  we can take the expectation on both sides of \eqref{eqn:CauchyBinet}.
  We get $1$ on the left-hand side and the sum of $\binom{d}{p}$ identical terms
  on the right-hand side.
  Hence, the average squared $p$-dimensional volume of the projection
  is $1 / \binom{d}{p}$, as claimed.
\end{proof}

We set $\DConstant{p}{d} = \moment{p}{d}{1} / \moment{p}{d}{2}$ and leave it
to the reader to verify that this agrees with \eqref{eqn:DConstant},
where $\DConstant{p}{d}$ is given in terms of Gamma functions as well as
double factorials.

\section{Tiling the Space of Point-Directions}
\label{sec:3}

We use the Delaunay mosaic to tile the space of \emph{point-direction pairs},
$\Rspace^d \times \Grass{p}{d}$.
Given a Delaunay mosaic of a set $A \subseteq \Rspace^d$, denoted $\Delaunay{A}$,
consider a $p$-dimensional cell, $\gamma \in \Delaunay{A}$,
and its dual $(d-p)$-dimensional Voronoi cell, $\gamma^* \in \Voronoi{A}$.
We define the \emph{$p$-tile} of $\gamma$ to consist of all
pairs $(x, L) \in \Rspace^d \times \Grass{p}{d}$ such that
$L+x$ has a non-empty intersection with $\gamma^*$,
and $x$ lies in the projection of $\gamma$ onto $L+x$:
\begin{align}
  \tile{\gamma}{\gamma^*}  &=  \{ (x,L) \in \Rspace^d \times \Grass{p}{d}
    \mid  x \in \proj{\gamma}{L+x}
          \mbox{\rm ~and~} (L+x) \cap \gamma^* \neq \emptyset \} .
\end{align}
The \emph{tiles} decompose the space $\Rspace^d \times \Grass{p}{d}$
in the sense that they cover the space of point-direction pairs
and their interiors are pairwise disjoint.
Since the detailed analysis of the boundaries is irrelevant for the
current work, we only prove a weaker statement.
\begin{lemma}[Uniqueness of Tile]
  \label{lem:uniqueness_of_tile}
  Let $A \subseteq \Rspace^d$ be locally finite with $\conv{A} = \Rspace^d$,
  and let $0 \leq p \leq d$.
  Then for almost every point-direction pair,
  $(x, L) \in \Rspace^d \times \Grass{p}{d}$,
  there exists a unique $p$-tile, $\tile{\gamma}{\gamma^*}$, that contains $(x, L)$.
\end{lemma}
\begin{proof}
  Take any point-direction pair, $(x, L)$.
  Assume without loss of generality that $x = 0$ is the origin
  and $L = \Rspace^p$ is a coordinate $p$-plane in $\Rspace^d$.
  Map each point $a \in A$ to the point $a' = \proj{a}{L} \in \Rspace^p$,
  and let $a'' = - \Edist{a}{a'}^2 \in \Rspace$ be its \emph{weight}.
  The weighted points define a weighted Voronoi tessellation
  and the corresponding weighted Delaunay mosaic; see e.g.\ \cite{Aur87b}.
  The mosaic is generically a simplicial complex and generally a polyhedral complex,
  which is geometrically realized in $\Rspace^p$ by drawing each cell, $\gamma$,
  as the convex hull of the points that generate the $p$-dimensional Voronoi cells
  sharing $\gamma^*$.
  Consider the cells in $\Voronoi{A}$ that have a non-empty intersection
  with $L$, write $\cScape{L}{A}$ for the collection of
  dual cells in $\Delaunay{A}$, and
  observe that $\cScape{L}{A}$ is the Voronoi scape of $L$.

  As proved in \cite{Sib80}, the weighted Voronoi tessellation is
  the intersection of $L$ with $\Voronoi{A}$ and, by duality,
  the weighted Delaunay mosaic is the orthogonal projection of $\cScape{L}{A}$ to $L$.
  If $L$ and $\Delaunay{A}$ are in general position,
  then all Delaunay cells in $\cScape{L}{A}$ project injectively to $L$,
  and the cells of dimension less than $p$ form a set of zero measure.
  If $\Delaunay{A}$ covers $\Rspace^d$,
  then the weighted Delaunay mosaic covers $\Rspace^p$.
  Hence, for almost all point-direction pairs, $(x, L)$,
  there is a unique Delaunay $p$-cell $\gamma$,
  such that $(x, L) \in \tile{\gamma}{\gamma^*}$, as claimed.
\end{proof}

The proof of the lemma gives some insight into the motivation for choosing
this particular tiling of the space of point-direction pairs.
We now compute the measure of a tile.
\begin{lemma}[Volume of Tile]
  \label{lem:volume_of_tile}
  The measure of $J = \tile{\gamma}{\gamma^*}$ is
  $\Vol{}{J}{} = \Vol{p}{\gamma}{} \, \Vol{d-p}{\gamma^*}{} / \binom{d}{p}$.
\end{lemma}
\begin{proof}
  The measure of the tile is the integral of $1$ over its pairs.
  Setting $x = y + z$, in which $y \in L$ and $z \in L^\bot$, the integral is
  \begin{align}
    \Vol{}{J}{}  &=  \int\nolimits_{L \in \Grass{p}{d}}
        \int\nolimits_{y \in L} \Indicator{y \in \proj{\gamma}{L}}
        \int\nolimits_{z \bot L} \Indicator{(L+z) \cap \gamma^* \neq \emptyset}
             \diff z \diff y \diff L
    \label{eqn:tilemusic1} \\
                 &=  \Vol{p}{\gamma}{} \, \Vol{d-p}{\gamma^*}{} \,
        \int\nolimits_{L \in \Grass{p}{d}} \cos^2 \Angle{L}{\gamma} \diff L ,
    \label{eqn:tilemusic2}
  \end{align}
  where we get \eqref{eqn:tilemusic2} by noticing that
  the innermost integral in \eqref{eqn:tilemusic1}
  is the $(d-p)$-dimensional volume of the projection of $\gamma^*$ to $L^\bot$,
  which is $t \, \Vol{d-p}{\gamma^*}{}$ with $t = \cos \Angle{L}{\gamma}$,
  and the middle integral is the $p$-dimensional volume of the projection
  of $\gamma$ to $L$, which is $t \, \Vol{p}{\gamma}{}$.
  Using \eqref{eqn:momentDef2} and \eqref{eqn:moment2}, we see that
  the integral in \eqref{eqn:tilemusic2} is $\moment{p}{d}{2}$.
\end{proof}

We take a closer look at the projection of a tile to $\Rspace^d$.
Let $(x, L)$ be a point-direction pair in $J = \tile{\gamma}{\gamma^*}$
with $\dime{\gamma} = p$.
There are points $u \in \gamma$ and $v \in \gamma^*$ such that
$x = \proj{u}{L+x}$ and $v = (L+x) \cap \gamma^*$.
Because of the right angle between the direction and the projection,
we have $\Edist{x}{u}^2 + \Edist{x}{v}^2 = \Edist{u}{v}^2$,
so $x$ lies on the smallest sphere that passes through $u$ and $v$.
Indeed, $u$ and $v$ define a $(d-1)$-dimensional set of
point-direction pairs, and the points of these pairs all lie on
the mentioned sphere.
Let $z_0 = z (\gamma, \gamma^*) = \affine{\gamma} \cap \affine{\gamma^*}$
and observe that the sphere defined by $u$ and $v$ also passes
through $z_0$.
Let $R_0 = R (\gamma, \gamma^*)$ be the maximum distance between a
point of $\gamma$ and a point of $\gamma^*$ and note that $R_0$ is
the radius of every largest empty sphere that passes through the
vertices of $\gamma$.
Since the diameter of the sphere spanned by $u$ and $v$ is
$\Edist{u}{v} \leq R_0$,
it follows that the ball with center $z_0$ and radius $R_0$ contains
this sphere.
This implies an upper bound on the size of the projection of the
tile to $\Rspace^d$, which we state formally for later reference.
\begin{lemma}[Projection of Tile]
  \label{lem:projection_of_tile}
  The projection of $J = \tile{\gamma}{\gamma^*}$ to $\Rspace^d$ is
  contained in the closed ball with center $z_0 = z(\gamma, \gamma^*)$
  and radius $R_0 = R(\gamma, \gamma^*)$.
\end{lemma}

It follows that the volume of the projection of the tile to $\Rspace^d$
is at most $R_0^d$ times the volume of the unit ball in $\Rspace^d$.
Since we assume the uniform probability measure on $\Grass{p}{d}$,
the same upper bound holds for the measure of the tile itself.

\section{Mixed Regularity}
\label{sec:4}

Taking the union of progressively more tiles,
we eventually cover all of $\Rspace^d \times \Grass{p}{d}$.
However, at each step some of the points miss some of the directions,
and which directions are covered depends on the mosaic.
For what follows, we require a mild regularity condition for this tiling.
For a set $\Omega \subseteq \Rspace^d$ we call a tile
a \emph{boundary tile of $\Omega$} if its projection to $\Rspace^d$
contains at least one point inside and at least one point outside $\Omega$.
\begin{definition}[Mixed Regularity]
  \label{def:mixed_regularity}
  Let $A \subseteq \Rspace^d$ be locally finite.
  We say that $A$ has the property of \emph{mixed regularity} if, for any $p$,
  the total measure of the boundary $p$-tiles of a $d$-ball of radius $R$ centered
  at the origin is $o(R^d)$.
\end{definition}

Note that $\conv{A} = \Rspace^d$ is necessary for $A$ to have
the mixed regularity property.
Indeed, if $\conv{A}$ does not cover $\Rspace^d$, then there exists an unbounded
Voronoi cell and thus a tile with infinite measure.
Motivated by the analysis in Section~\ref{sec:3},
we give some sufficient conditions for a set $A \subseteq \Rspace^d$
to satisfy the mixed regularity property.
\begin{lemma}[Sufficient Conditions]
  \label{lem:sufficient_conditions}
  A locally finite set $A \subseteq \Rspace^d$ has the mixed regularity property
  if one of the following holds:
  \begin{enumerate}
    \item The radii of all circumspheres of top-dimensional Delaunay cells are bounded.
    \item Each ball in $\Rspace^d$ of radius greater than $R_0$ contains a point of $A$.
    \item There is a function $g(R) = o(R)$ such that every ball
      of radius $g(R)$ that intersects the $d$-ball of radius $R$
      centered at the origin contains at least one point of $A$.
  \end{enumerate}
\end{lemma}

Conditions 1 and 2 are equivalent, while the last one is weaker.
We finish the section with an application to Poisson point processes.
\begin{lemma}[Mixed Regularity in Expectation]
  \label{lem:mixed_regularity_in_expectation}
  A stationary Poisson point process $A \subseteq \Rspace^d$
  satisfies the mixed regularity property in expectation; that is:
  the total \emph{expected} measure of the boundary tiles of a $d$-ball
  of radius $R$ centered at the origin is $o(R^d)$.
\end{lemma}

\begin{proof}
  For each boundary tile, its Delaunay cell (which is almost surely a simplex)
  is a face of a top-dimensional Delaunay simplex
  whose circumsphere intersects the boundary of the ball.
  In \cite[Appendix A]{ENR17} it is proved that the total number
  of such spheres is $o(R^d)$.
  A minor modification of that proof suffices to show that the total volume
  of these balls is $o(R^d)$.
  Together with Lemma~\ref{lem:projection_of_tile}, this implies the
  claimed bound.
\end{proof}

\section{Mixed Volume}
\label{sec:5}

Call $\Vol{p}{\gamma} \, \Vol{d-p}{\gamma^*}$ the \emph{mixed volume} of a $p$-cell
$\gamma \in \Delaunay{A}$ and its dual $(d-p)$-cell $\gamma^* \in \Voronoi{A}$.
We note that this concept is related to a particular decomposition of $\Rspace^d$,
as we now explain.
Given $A \subseteq \Rspace^d$, the $d$-dimensional cells of the \emph{mixed complex}
defined in \cite{Ede99} are translates of the products
$\frac{1}{2} \gamma \times \frac{1}{2} \gamma^*$.
The $d$-dimensional volume of this cell is
$\Vol{d}{\frac{1}{2} \gamma \times \frac{1}{2} \gamma^*}{}
    = \Vol{p}{\gamma}{} \, \Vol{d-p}{\gamma^*}{} / 2^d$.
As proved in \cite{Ede99}, the cells in the mixed complex have
pairwise disjoint interiors and they cover $\Rspace^d$.
Assuming the mixed regularity property, this implies that, up to a lower order term,
the cells for $p = 0$ cover a fraction of $1 / 2^d$ of the ball of radius $R$
centered at the origin.
By symmetry, this is also true for $p = d$.
We continue with a generalization of these bounds to dimension $p$ between $0$ and $d$.
\begin{corollary}[Mixed Volumes]
  \label{cor:mixed_volumes}
  Let $A \subseteq \Rspace^d$ have the mixed regularity property.
  For any $0 \leq p \leq d$, the sum of the mixed volumes
  over all $p$-dimensional Delaunay cells contained in a ball of radius $R$ is
  $\nu_d \tbinom{d}{p} R^d + o(R^d)$,
  in which $\nu_d$ is the volume of the unit ball in $\Rspace^d$.
\end{corollary}
\begin{proof}
  Let $B(R)$ be the ball with radius $R$ centered at the origin.
  Set $\Bcal{p}{R} = B(R) \times \Grass{p}{d}$,
  let $M_p (R)$ be the smallest collection of $p$-tiles whose union
  contains $\Bcal{p}{R}$, and let $\partial M_p(R)$ be the boundary tiles of $B(R)$.
  Clearly,
  \begin{align}
    M_p(R) \setminus \partial M_p(R)  &\subseteq  \Bcal{p}{R}  \subseteq  M_p (R) .
  \end{align}
  We note, that if a tile, $J = J(\gamma, \gamma^*)$, contains a point inside the ball,
  then either $\gamma$ is inside the ball, or $J$ is a boundary tile.
  Indeed, for every point $x\in \gamma \setminus B(R)$, there is a direction $L$,
  such that $L+x$ intersects $\gamma^*$, hence $(x,L) \in J(\gamma, \gamma^*)$.
  By Lemma~\ref{lem:volume_of_tile}, the measure of this tile
  is $\Vol{p}{\gamma} \Vol{d-p}{\gamma^*} / \tbinom{d}{p}$ and,
  by the mixed regularity property,
  the measures of the tiles corresponding to Delaunay cells inside the ball
  sum up to $\Vol{d}{B(R)}(1 + o(1))$.
  Multiplying by $\tbinom{d}{p}$ concludes the proof.
\end{proof}

\section{Average and Expected Distortion}
\label{sec:6}

Recall that the Voronoi scape of a set $\Omega \subseteq \Rspace^d$ and
a locally finite set $A \subseteq \Rspace^d$ consists of all Delaunay cells
whose dual Voronoi cells have a non-empty intersection with $\Omega$.
Assuming $\Omega$ is $p$-dimensional, only intersections with
$(d-p)$-dimensional Voronoi cells are relevant.
Generically, such an intersection has a finite \emph{multiplicity},
and we define the \emph{Voronoi scape}, denoted $\cScape{\Omega}{A}$,
as the multiset of $p$-cells, $\gamma \in \Delaunay{A}$,
in which $\gamma$ appears the multiplicity
of the intersection between $\Omega$ and $\gamma^*$ times.
We are ready to prove the main result of the paper.
\begin{theorem}[Average Volume]
  \label{thm:average_volume}
  Let $A \subseteq \Rspace^d$ have the mixed regularity property,
  and let $\Omega$ be a $p$-dimensional rectifiable set in $\Rspace^d$.
  The average volume of $\cScape{\Omega}{A}$ over all congruent copies
  of $\Omega$ inside the $d$-ball with radius $R$ centered at the origin is
  $\Vol{p}{\Omega} (\DConstant{p}{d} + o(1))$ as $R$ goes to infinity.
\end{theorem}
\begin{proof}
  We start with the Crofton formula \cite[Formula (5.7)]{ScWe08},
  which states that
  \begin{align}
    \int_{Q \in \Grass{d-p}{d}} \int_{y \in Q^\bot}
      \chi((Q + y) \cap \Omega) \diff y \diff Q
    &= \moment{p}{d}{1} \Vol{p}{\Omega},
  \end{align}
  in which $\chi((Q+y) \cap \Omega)$ is the almost always finite
  multiplicity of the intersection between $Q+y$ and $\Omega$.
  Fixing a $(d-p)$-dimensional moving polyhedron, $P$, we notice
  that for each intersection point of $Q$ and $\Omega$,
  the total measure of congruent copies of $P$
  that intersect $\Omega$ at this point is $\Vol{d-p}{P}$.
  Hence, the total measure of intersection points over all congruent copies
  of $P$ is
  \begin{align}
    \int_{P' \cong P} \chi(P' \cap \Omega) \diff P'
      &= \Vol{d-p}{P} \int_{Q \in \Grass{d-p}{d}}
         \int_{y \in Q^\bot} \chi((Q+y) \cap \Omega) \diff y \diff Q \\
      &= \Vol{d-p}{P} \Vol{p}{\Omega} \moment{p}{d}{1}.
  \end{align}
  Taking $P = \gamma^*$ and moving $\Omega$ instead of the polyhedron,
  we see that the total measure of intersection points of congruent copies
  of $\Omega$ with $\gamma^*$ is $\Vol{d-p}{\gamma^*} \Vol{p}{\Omega} \moment{p}{d}{1}$.

  A $p$-cell $\gamma \in \Delaunay{A}$ belongs to the Voronoi scape
  of a congruent copy $\Omega'$ of $\Omega$ precisely $\chi(\Omega' \cap \gamma^*)$ times,
  and we just computed the total value of this quantity.
  The total contribution of $\gamma$ to the $p$-volume of the Voronoi scapes
  of the congruent copies of $\Omega$ is thus
  $\Vol{p}{\gamma} \Vol{d-p}{\gamma^*} \Vol{p}{\Omega} \moment{p}{d}{1}$.
  We get the final result by dividing the total contribution of the Delaunay
  cells inside the ball---which we compute using
  Corollary \ref{cor:mixed_volumes}---by the total measure
  of the congruent copies inside the ball:
  \begin{align}
    \frac{ \Vol{d}{B(R)} \Vol{p}{\Omega} \moment{p}{d}{1} / \moment{p}{d}{2} }
         { \Vol{d}{B(R)}(1 + o(1)) }
      &=  \Vol{p}{\Omega} (\DConstant{p}{d} + o(1)) ,
  \end{align}
  in which we use \eqref{eqn:moment1} and \eqref{eqn:moment2} to get the
  expression on the right.
\end{proof}

We finish with stating the answer to the original question that motivated
the work reported in this paper.
We showed in Section~\ref{sec:4} that the stationary Poisson point process
has the mixed regularity property in expectation,
which allows us to repeat all results while
adding the expectation to all quantities.
By the isometry invariance of the process,
for any set $\Omega$, the expected volume of $\cScape{\Omega}{A}$
does not depend on the position of $\Omega$.
Exchanging the expectation and the average inside the ball of radius $R$
centered at the origin and letting $R$ go to infinity,
we arrive at probabilistic versions of
Theorem~\ref{thm:average_volume}.
\begin{theorem}[Expected Volume]
  \label{thm:expected_volume}
  Let $A \subseteq \Rspace^d$ be a stationary Poisson point process
  with intensity $\rho > 0$, and let $\Omega$ be a compact
  rectifyable $p$-manifold in $\Rspace^d$.
  Then the expected $p$-dimensional volume of the Voronoi scape
  of $\Omega$ is $\DConstant{p}{d} \Vol{p}{\Omega}$.
\end{theorem}

\section{Discussion}
\label{sec:8}

The main contribution of this paper is a complete analysis
of the average and expected distortion of $p$-dimensional Voronoi scapes
in $\Rspace^d$, for $0 \leq p \leq d$.
For $p = 1$, these scapes are known as Voronoi paths,
for which the expected distortion has been studied but was known
only in $\Rspace^2$; see \cite{BTZ00}.
A useful insight from our analysis is that the expected distortion
for a stationary Poisson point process is the average distortion
of a general locally finite point set.
We make crucial use of this insight in the proof of our results.
\emph{Can these results be extended to other measures,
such as notions of curvature, for example?}
The proof of Theorem~\ref{thm:average_volume}
suggests that this extension
would require a detailed analysis of the Crofton formula.
Insights in this direction could be helpful in using the Voronoi scape
to measure otherwise difficult to measure shapes \cite{BEEN20}.

In our analysis, the properties that make a mosaic a Delaunay mosaic
are not used other than in the quantification of the mixed regularity
property for locally finite sets.
Indeed, we only need a pair of dual complexes in which dual cells
are orthogonal to each other,
a property that holds also for the generalizations of Voronoi
tessellations and Delaunay mosaics to points with real weights;
see e.g.\ \cite{Aur87b}.

\subsection*{Acknowledgements}

{\footnotesize
  The authors thank Ranita Biswas and Tatiana Ezubova for the collaboration
  on an experimental project that motivated the work reported in this paper.
}


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

Scattered radiation conveys special information regarding both the emission and scattering regions.
In particular, linear polarization may develop in an anisotropic scattering geometry and the relative motion
between the emission and scattering regions gives rise to complicated line profiles including multiple-peak
structures and broad wings. An excellent example is provided by spectropolarimetric observations
of Seyfert 2 galaxies exhibiting broad lines in the linearly polarized fluxes, lending strong support
to the unification model of active galactic nuclei
\citep{miller90, tran10}.

Interaction of electromagnetic radiation with an atomic electron may be classified into Rayleigh
and Raman scattering \citep[e.g.][]{sakurai67}. Rayleigh scattering refers to the elastic process where the scattered photon
has the same wavelength as the incident one. Otherwise, we have Raman scattering in which the initial and final electronic states differ so that the scattered photon emerges with the energy difference, which is enforced by the law of
energy conservation. Raman spectroscopy is particularly useful in revealing the complicated energy level structures
of a molecule in chemistry.

In astrophysics,
\cite{nussbaumer89} provided a pioneering discussion on Raman scattering with atomic hydrogen, introducing
a new spectroscopic tool to diagnose gaseous emission nebulae including symbiotic stars and active galactic
nuclei. Relevant Raman scattering processes start with a far UV photon more energetic than Ly$\alpha$
incident on an hydrogen atom in the ground state and end with an outgoing photon with energy less than
that of Ly$\alpha$ leaving behind the hydrogen atom in the $2s$ state.
They introduced basic atomic physics of Rayleigh and Raman scattering illustrating the cross sections
and presented a number of candidate far UV spectral lines that may result in detectable Raman-scattered features
\citep[e.g.,][]{saslow69}.
One notable point is that the branching ratio of Raman and Rayleigh scattering of far UV radiation near Ly$\beta$
is approximately 0.15 so that most Ly$\beta$ photons are Rayleigh (or resonantly) scattered several times before
they are converted into H$\alpha$ photons to escape from the thick neutral region.

\cite{schmid89} identified the broad emission features at 6825 \AA~and 7082 \AA~ found
in about half of symbiotic stars as Raman scattered features of  O VI $1032$ and $1038$ \AA~ emission
lines, respectively \citep{akras19}. Symbiotic stars are wide binary systems composed of a hot white dwarf and
a mass-losing red giant. Hydrodynamical studies suggest that
some fraction of slow stellar wind is gravitationally captured to form an accretion disk \citep{devalborro17,chen17,saladino18}, where
O VI$\lambda\lambda$1032 and 1038 are important coolants.

Considering very small Raman scattering cross sections $\sim 10^{-23}
{\rm\ cm^2}$ for O VI$\lambda\lambda$1032 and 1038 doublet lines, the operation of Raman scattering
requires a special condition that a very thick neutral
region is present in the vicinity of a strong O VI emission region. This special condition is ideally met in symbiotic
stars, where a thick neutral region surrounding the giant component is illuminated by strong far UV radiation
originating from the nebular region that may be identified with the accretion flow onto the hot component.

Raman-scattered O VI features in symbiotic stars exhibit
complicated profiles with multiple peaks separated by $\sim 30-50{\rm\ km\ s^{-1}}$ indicative
of the O VI emission regions with physical dimension of $\sim 1{\rm\ au}$ \citep[e.g.][]{shore10,heo15, lee19}.
Because Raman and Rayleigh scattering sufficiently off resonance shares the same scattering phase function
as Thomson scattering\citep{schmid95,yoo02,chang17}, strong linear polarization may develop in an anisotropic scattering geometry.
\cite{harries96} conducted spectropolarimetric observations of many symbiotic stars to show that Raman-scattered
O VI features are strongly polarized. They also found that most Raman-scattered O VI features show polarization
flip in the red wing part, where the polarization develops nearly perpendicularly to the direction along which
the main part is polarized. \cite{lee99} proposed that the polarization flip is closely associated with the bipolar
structure of symbiotic stars \citep[see also][]{heo16}.

Raman scattering plays an interesting role of redistributing far UV radiation near Ly$\beta$ and Ly$\gamma$
into near H$\alpha$ and H$\beta$, respectively. He~II being a single electron ion with $Z=2$, the transition
to $n=2$ from an energy level from $n=2k, k>1$ gives rise to emission lines with wavelengths slightly shorter than
those of H I Lyman series $k \to 1$, for which the cross sections for Rayleigh and Raman scattering
are conspicuously large. \cite{vangroningen93} found Raman-scattered He~II features near H$\beta$ in
the symbiotic nova RR~Telescopii. \cite{pequignot97} discovered the same spectral feature
in the young planetary nebula NGC~7027, which constitutes the first discovery of a spectral feature formed through Raman
scattering with atomic hydrogen in planetary nebulae.
Subsequently, Raman-scattered He~II at $6545$ \AA~ has been detected in several symbiotic stars \citep{birriel04,jung04,sekeras15}
and in the young planetary nebulae NGC~6302, IC~5117, NGC~6790, NGC~6886, and NGC~6881 \citep{groves02, lee01,lee06,lee09,letter20}.

It is particularly notable that the
case B recombination theory allows one to deduce the strengths of incident far UV He~II lines \citep{storey95}, yielding faithful
estimates of Raman conversion efficiencies. This is extremely useful in the measurement of
H~I content in symbiotic stars and young planetary nebulae, which makes Raman spectroscopy a totally new
approach to probing the mass loss processes occurring in the late stage of stellar evolution
\citep{lee06,choi20}.
Furthermore, H$\alpha$ and H$\beta$ in symbiotic stars and young planetary nebulae often display fairly
extended wings that may indicate the
presence of fast tenuous stellar wind \citep[e.g.,][]{arrieta02}. Broad wings around Balmer lines may also arise via Raman scattering far-UV
continuum around Lyman lines, which requires further investigation \citep{lee00a,yoo02,chang18}.

Additional examples include Raman scattering of C II$\lambda\lambda$1036 and 1037 forming
optical features at 7023 \AA\ and 7054 \AA, which
were reported in the symbiotic nova V1016~Cygni by \cite{schild96}. \cite{dopita16} investigated
the H II regions in the Orion Nebula (M42) and five H II regions in the Large and Small
Magellanic Clouds to discover Raman-scattered features at 6565 \AA\ and 6480 \AA\ formed through Raman scattering of O I $\lambda$1025.76 and Si II $\lambda$1023.70, respectively.

The Monte Carlo approach is an efficient numerical technique to describe radiative transfer in various regions having
	 a dust component\citep[e.g.][]{celnikier74,seon15} and a molecular component\citep[e.g.][]{brinch10}.
	 A similar approach has been applied to radiative transfer of an
	electron scattering \citep[e.g.][]{angel69,seon94}. Ly$\alpha$ deserves special attention being
	characterized by large scattering optical depth \citep[e.g.][]{eide18,seon20}.

\cite{schmid92} performed  Monte-Carlo simulations to investigate the formation of Raman-scattered O VI
features in an expanding H~I region with an assumption that a given line photon has an invariant scattering cross
section as it propagates through the H~I medium. Similar studies were presented by \cite{lee97b}, who
adopted a density matrix formalism to determine the physical information of scattered radiation including polarization.
\cite{chang15} investigated the formation of broad Balmer wings near H$\alpha$ and H$\beta$
in the unification scheme of active galactic nuclei and presented quantitatively
the asymmetry of the wings formed in neutral regions with extremely high H I column density $\sim 10^{23} \rm\ cm^{-2}$.

In this paper, we introduce a new grid-based Monte Carlo code entitled "Sejong Radiative Transfer for Rayleigh
and Raman Scattering ({\it STaRS}),
We also present our test of the code by revisiting the formation of Balmer wings and Raman O VI features
through Raman scattering with atomic hydrogen.

\section{Grid Based Monte Carlo Simulation}\label{sec:method}

In this section, we describe 'STaRS' and discuss the basic atomic physics of Rayleigh and Raman scattering with atomic hydrogen.
Fig.~\ref{fig:scattering} is a schematic illustration of a few
representative transitions pertaining to the two types of scattering. Thus far, detected Raman-scattered
features are limited to those associated with the final de-excitation to $2s$ state.
The second order time-dependent perturbation theory is used to compute the scattering cross sections
known as the Kramers-Heisenberg formula \citep[e.g.,][]{bethe67, sakurai67, saslow69}.

In our grid based Monte Carlo code, we divide
the region of interest into a large number of small cubes or cells in the Cartesian coordinate system, where each cell is characterized by uniform physical
properties in the three dimensional space.
Here, the uniform physical properties include
H I number density $n_{{\rm HI},G}$ and the velocity ${\bf v}_{G}$.
We  assign the emissivity $j_{e}$ to each cell and
generate an initial photon using $j_{e}(\lambda_i,x,y,z)$, where $\lambda_i$ is the wavelength of the initial photon.
STaRS is mainly written in {\it Fortran} with Message Passing Interface implemented for parallel computing.
We also adopt the shared memory technique by {\it intel MPI}.

Fig.~\ref{fig:flow_chart} shows a flow chart for radiative transfer simulations using STaRS.
A simulation starts with the setup of the scattering geometry by assigning to each cell appropriate
physical conditions. Initial far UV photons are generated in accordance with our prescription of $j_{e}$.
Each photon is tagged with the information of
the unit wavevector $\bf{\hat k}$, the position vector $\bf r$, the wavelength $\lambda$,
and the density matrix $\rho$ composed of the Stokes parameters, $I$, $Q$, $U$, and $V$.
The free path $d$ for next scattering position $\bf r'$ is computed by transforming the scattering optical
depth into the physical depth.

Decision is made whether the photon escapes from the region or is scattered into another direction.
If the next scattering position $\bf r'$ is outside the scattering region, the photon is assumed to
reach the observer as a far UV photon. Otherwise, we generate a new photon using the scattering
phase function for an electric dipole process and determine the scattering type.
We assume that the region is transparent to Raman-scattered photons. The scattered photon escapes
from the scattering region if the scattering is Raman. Otherwise, the Rayleigh scattered photon
is regarded as an incident photon propagating to a new scattering position.
The procedure is repeated until escape.

In Sec.~\ref{sec:geometry} and \ref{sec:source}, we provide
more detailed descriptions on the scattering geometry and generation of initial far UV photon prescribed
by $j_{e}$. In Sec.~\ref{sec:path}, we describe the computation of a free path $d$.
We discuss the basic properties of the scattered photons in Sec.~\ref{sec:scattering}.

\begin{figure}
	\includegraphics[width=90mm]{scattering.png}
	\caption{Schematic illustrations of energy levels and electronic transitions associated with Raman and Rayleigh scattering with a hydrogen atom of far UV electromagnetic radiation.
	}
	\label{fig:scattering}
\end{figure}

\begin{figure}
	\includegraphics[width=90mm]{flow_chart.png}
	\caption{A flow chart of STaRS.
	}
	\label{fig:flow_chart}
\end{figure}

\subsection{Geometry : Scattering Region}\label{sec:geometry}

The medium for radiative transfer through Raman and Rayleigh scattering
corresponds to a thick H I region that is easily found in the slow stellar wind from a red giant \citep{lee97b,lee19}.
No consideration on the thermal motion of neutral hydrogen is given
because the variation of the cross section and the branching ratio in the scale of thermal speed
is negligible \citep{chang15,chang18}.

In our simulation, we divide the scattering region into a large number of cells having a fixed size.
Each cell is identified by a set of three indices $(iX, iY, iZ)$ with $iX, iY$ and $iZ$ running from 1 to
$N_X, N_Y$ and $N_Z$.
If we denote by $X_{min}$ and $X_{max}$ the range of $x$ coordinates for the scattering region,
the boundary of the $iX$th cell on the $x$ axis is given by
\begin{equation}
X_G(iX) = X_{min} + \frac{X_{max} - X_{min}}{N_X} (iX - 1).
\end{equation}
We also define $Y_G(iY)$ and $Z_G(iZ)$ in a similar way.

Therefore, for any point $P(x,y,z)$ in the cell identified with $(iX,iY,iZ)$, we have
\begin{eqnarray}
X_G(iX) < x < X_G(iX+1) \\
\nonumber
Y_G(iY) < y < Y_G(iY+1) \\
\nonumber
Z_G(iZ) < z < Z_G(iZ+1).
\end{eqnarray}
Each cell is also identified by its center point, whose coordinates
are given by the following relations
\begin{eqnarray}
X_C(iX) = \frac{X_G(iX) + X_G(iX+1)}{2} \\
\nonumber
Y_C(iY) =  \frac{Y_G(iY) + Y_G(iY+1)}{2} \\
\nonumber
Z_C(iZ) =  \frac{Z_G(iZ) + Z_G(iZ+1)}{2}.
\end{eqnarray}

\subsection{Emission Source : Initial Photon}\label{sec:source}

The operation of Raman scattering with atomic hydrogen requires the coexistence of a strong far UV source
and a thick neutral region. In the case of symbiotic stars and young planetary nebulae, a strong far UV
emission region is formed near an accreting white dwarf or a hot central star and a thick neutral region
is also present in association with the mass loss of a giant star. In the simulation,
the emissivity $j_{e}(\lambda,x,y,z)$ is prescribed as a function of wavelength and position.
Regarding $j_{e}$ as the normalized probability density function,
we pick a wavelength $\lambda_i$ and a starting
position $(x_i,y_i,z_i)$ of an initial photon in accordance with $j_{e}(\lambda,x,y,z)$.

We find the spatial index $(iX,iY,iZ)$ from $(x_i,y_i,z_i)$ to determine the cell that contains the starting
position. For simplicity, we assume that the initial photon is completely unpolarized and that the unit
wavevector of the initial photon is selected from an isotropic distribution.
Thus, the initial unit wavevector ${\bf \hat{k}}=(k_x,k_y,k_z)$ is obtained using two uniform random
numbers $r_1$ and $r_2$ between 0 and 1 with the following prescription
\begin{eqnarray}
\mu &=& \cos\theta = 2r_1 - 1 \\
\nonumber
\phi &=& 2\pi r_2 \\
\nonumber
k_x &=& \sin\theta \cos \phi \\
\nonumber
k_y &=& \sin\theta \sin \phi \\
\nonumber
k_z &=& \cos\theta.
\end{eqnarray}\label{eq:wavevector}
Because scattered radiation is polarized in an anisotropic geometry, it is important to carry the polarization information.
The four Stokes parameters $(I,Q,U,V)$ are necessary to describe the polarization state. An equivalent way is provided by considering the $2\times2$ density matrix defined by
\begin{equation}
\rho =
\begin{bmatrix}
{(I + Q)/2} & {(U + iV)/2} \\
{(U - iV)/2} & {(I - Q)/2} \\
\end{bmatrix}
.
\end{equation}.

In our simulation, the Stokes parameter $V$, representing circular polarization, is always set to zero because
no circular polarization develops from initially unpolarizd photons in electric dipole processes associated with Raman and Rayleigh scattering.
Initially unpolarized photons are described by a simple $\rho$ is given by
\begin{eqnarray}
\rho_{11} &=& 0.5 \\
\nonumber
\rho_{22} &=& 0.5 \\
\nonumber
\rho_{12} &=& \rho_{21} = 0\ .
\end{eqnarray}

The wavelength measured by an observer lying on the photon path is kept in the simulation. We assume that
the emitters are in random motion and also subject to the bulk motion associated with the grid.
If we let the velocity of the emitter be ${\bf v}_{\rm emit}$
including bulk and random velocities,
the wavelength in the grid frame $\lambda_g$ is given by
\begin{equation}
\lambda_g = \lambda_i \left(  1 - {\frac{{\bf v}_{\rm emit} \cdot {\bf\hat k}}{c} } \right)
  \left(  1 + {\frac{{\bf v}_{\rm G} \cdot {\bf\hat k}}{c}} \right
),
\end{equation}
where ${\bf v}_{\rm G}$ is the velocity of the cell $(iX,iY,iZ)$.

\subsection{Journey of Photon: Optical Depth and Free Path }\label{sec:path}

\begin{figure*}
	\includegraphics[width=180mm]{cross.eps}
	\caption{The total scattering cross section(black),
	the branching ratio to $n=2$ (blue) and $n=3$ (red)
	computed in \cite{chang15}.
	}
	\label{fig:cross}
\end{figure*}

The generation of an initial photon is followed by the estimate of the free optical depth $\tau$ given by
\begin{equation}
\tau = - \ln r,
\end{equation}
where $r$ is a uniform random number between $0$ and $1$.
We adopt the method in \cite{seon09} to compute the free path and the next scattering position in grid-based
geometry. In this section, we provide a brief description to compute the free path and determine the scattering position.

If the starting position ${\bf r}$ of the photon is in cell $A$, we measure the distance to the boundary
of the cell $d_A$ along the photon ray from the starting position. With $d_A$ we define the
scattering optical depth $\tau_A$ to the cell boundary by
\begin{equation}\label{eq:tau1}
\tau_A = \sigma_{\rm tot}(\lambda_A) n_{{\rm HI},A} d_A,
\end{equation}
where $\lambda_A$ is the wavelength in cell $A$.
Fig.~\ref{fig:cross} shows
the total scattering cross section $\sigma_{\rm tot}(\lambda)$ as a function of wavelength.
If $\tau_A>\tau$, then
the next scattering position ${\bf r}'$ is found in cell $A$ as follows
\begin{equation}
{\bf r'} = {\bf r} + \left( \frac{\tau}{\tau_A} \right) d_A {\bf \hat k}  .
\end{equation}

In the opposite case where $\tau_A < \tau$, the photon enters the neighboring cell $B$.
In this case, the same problem is obtained if we regard the entry point
${\bf r}_e={\bf r}+d_A{\bf\hat k}$
as the new starting point of the photon with a new free optical path $\tau'$ reduced by $\tau_A$, or
\begin{equation}
{\tau}' = {\tau} - {\tau_A}.
\end{equation}
It should be noted that we are dealing with radiative transfer in a medium in motion. Therefore,
cell $B$ may move with a velocity different from that of cell $A$, in which case the photon wavelength along its
propagation direction may change on entering cell $B$ from cell $A$. Denoting by ${\bf v}_{\rm G,A}$ and ${\bf v}_{\rm G,B}$
the velocities of cells $A$ and $B$, respectively, we have
\begin{equation}
\lambda_{B} = \lambda_A \left( 1 - {\frac{{\bf v}_{\rm G,A} \cdot {\bf\hat k}}{c}} \right)
\left( 1 + \frac{{\bf v}_{\rm G,B} \cdot {\bf\hat k}}{c} \right).
\end{equation}

Iterations from Eq.~\ref{eq:tau1} are made with necessary updates ${\bf r}={\bf r}_e$, $\tau=\tau'$, and $\lambda_A = \lambda_B$ and
new naming of cell $B$ as cell $A$ until we have $\tau_A>\tau$.
Fig.~\ref{fig:grid} shows a schematic illustration of this procedure.
In cases when a neighboring cell $B$ may not exist and $\bf r$ is outside the geometry,
the photon is regarded as Rayleigh-escaped.

\begin{figure}
	\includegraphics[width=90mm]{grid.png}
	\caption{Schematic illustration of computing the next scattering position in the grid-based simulation.
	}
	\label{fig:grid}
\end{figure}

\subsection{Rayleigh and Raman Scattering}\label{sec:scattering}

Raman scattering with atomic hydrogen with the initial and final electronic states being $1s$ and $2s$
shares the same the scattering phase function with Rayleigh scattering with atomic hydrogen in $1s$
\citep[e.g.][]{chang15}. It is also notable that the scattering is sufficiently far from resonance,
the scattering phase function is the same as that of Thomson scattering.
In this section, we describe the density matrix formalism with which
the polarization and the unit wavevector $\bf \hat k'$ of the scattered photon is determined  \citep{ahn15,chang17}.
We also describe the wavelength conversion and the line broadening associated with Raman scattering.

According to the density matrix formalism, the probability density of
the scattered wavevector ${\bf\hat k'}=(\sin\theta'\cos\phi',\sin\theta'\sin\phi', \cos\theta')$
is given by
\begin{equation}
I'({\bf\hat k'}) = \rho'_{11} + \rho'_{22},
\end{equation}
where $\rho'_{ij}$ is defined by
\begin{equation}
\rho'_{ij}=\sum_{kl=1,2}(\hat{\epsilon}'_i\cdot \hat{\epsilon_k})\rho_{kl}(\hat{\epsilon_l}\cdot \hat{\epsilon_l}').
\end{equation}
Here, $\hat{\epsilon}_{1,2}$ and $\hat{\epsilon}_{1,2}'$ are
the polarization basis vector associated with  $\bf \hat{k}$ and $\bf \hat{k}'$, respectively.
Specifically, $\hat\epsilon_1=(-\sin\phi,\cos\phi,0)$ and
$\hat\epsilon_2=(\cos\theta\cos\phi,\cos\theta\sin\phi, -\sin\theta)$ so that $\hat\epsilon_1$ represents
polarization in the direction perpendicular to the plane spanned by the photon wavevector and the
$z-$axis.

The components of the density matrix associated with the scattered radiation are related to those
of incident radiation by
\begin{eqnarray}
\rho'_{11} &=& (\cos^2\Delta\phi)\rho_{11}
\nonumber\\
&-& (\cos\theta \sin2\Delta\phi)\rho_{12}
\nonumber\\
&+& (\sin^2\Delta\phi \cos^2\theta)\rho_{22}
\nonumber\\
\rho'_{12} &=& (\frac{1}{2} \cos\theta'\sin2\Delta\phi)\rho_{11}
\nonumber\\
&+& (\cos\theta \cos\theta' \cos2\Delta\phi + \sin\theta \sin\theta' \cos\Delta\phi)\rho_{12}
\nonumber\\
&-& \cos\theta(\sin\theta \sin\theta' \sin\Delta\phi + \frac{1}{2} \cos\theta \cos\theta' \sin2\Delta\phi)\rho_{22}
\nonumber\\
\rho'_{22} &=& (\cos^2\theta' \sin^2\Delta\phi)\rho_{11}
\nonumber\\
&+& \cos\theta'(2\sin\theta \sin\theta' \sin\Delta\phi + \cos\theta \cos\theta' \sin2\Delta\phi)\rho_{12}
\nonumber\\
&+& (\cos\theta \cos\theta' \cos\Delta\phi + \sin\theta \sin\theta')^2 \rho_{22},
\end{eqnarray}
where $\Delta \phi = \phi'-\phi$.
In the code, a selection is made for ${\bf\hat k}'$ from an isotropic distribution to compute $I'$. A new
random deviate $r'$ is compared with $I'$. If $r'<I'$, then the selection of ${\bf\hat k}'$ is accepted.
Otherwise, the process is iterated until the acceptance is obtained.
The scattered photon now becomes a new incident photon propagating to a new scattering site. This is
followed by necessary updates for the wavevector and the density matrix given by
\begin{eqnarray}
{\bf \hat k} &=& {\bf \hat k'} \\
\nonumber
\rho_{ij} &=& \rho'_{ij}.
\end{eqnarray}

The scattering type is determined by the branching ratio.
If an incident far UV photon is more energetic than Ly$\gamma$, the final states
available for Raman scattering include $2s$, $3s$ and $3d$ states.
For photons near Ly$\gamma$, the total cross section is given by
\begin{equation}
\sigma_{\rm tot} = \sigma_{1s} + \sigma_{2s} + \sigma_{3s+3d},
\end{equation}
where $\sigma_{1s}$, $\sigma_{2s}$, and $\sigma_{3s+3d}$ are
the scattering cross sections corresponding to the final states $1s$, $2s$, and $3s+3d$, respectively.
Rayleigh branching ratio $BR1$ is given by
\begin{equation}
BR1 = \frac{\sigma_{1s}}{\sigma_{\rm tot}}.
\end{equation}
In a similar way,
Raman branching ratios corresponding to the final energy levels $n=2$ and $3$ are
\begin{equation}
BR2 = \frac{\sigma_{2s}}{\sigma_{\rm tot}}, \quad BR3 = \frac{\sigma_{3s+3d}}{\sigma_{\rm tot}},.
\end{equation}
In Fig.~\ref{fig:cross}, the blue and red solid lines represent $BR2$ and $BR3$, respectively.

The energy difference between the incident and Raman-scattered photons is the same as that
between the initial and final atomic states. This is translated into the relation between
the wavelengths $\lambda$ and $\lambda'$ of the Raman-scattered and incident photons, respectively,
which is given by
\begin{equation}\label{eq:wavelength}
\frac{1}{\lambda} = \frac{1}{\lambda'} + \frac{1}{\lambda_{\rm res}},
\end{equation}
where $\lambda_{\rm res}$ is the wavelength corresponding to the energy difference between
the initial and final states.

One very important aspect in Raman scattering can be found in the conspicuous change in line width.
Differentiating Eq.~(\ref{eq:wavelength}), we have
\begin{equation}
\frac{d\lambda'}{\lambda'}= \left( \frac{\lambda'}{\lambda} \right)
\left( \frac{d\lambda}{\lambda} \right),
\end{equation}
from which we see immediately that the line width of Raman-scattered feature is
broadened by the factor $(\lambda'/\lambda)^2$.
For example, a typically observed line width of Raman-scattered O VI at 6825 \AA amounts to $\sim 30$~\AA\
whereas the far UV parent line O VI$\lambda\lambda$1032 exhibits a line width $\sim 1$~\AA\ in many symbiotic
stars.

Due to the line broadening effect, far UV radiation near Lyman series of hydrogen will be considerably diluted
and redistributed around Balmer emission lines to appear as broad wings \citep[e.g.][]{yoo02,chang15,chang18}.
Another important consequence of the line broadening effect is that the line profiles of Raman-scattered
features mainly reflect the relative motion between the far UV emitters and the neutral scatterers and quite
independent of the observer's line of sight \citep{heo16,choi20}.

\section{Code Test}

As a check of our code, we present our simulation results for two exemplary cases.
The first example is a static spherical H I region surrounding a far UV
continuum source located at the center, in which Balmer wings are formed through Raman scattering.
Analytic solutions are available for this case,
against which our result obtained from STaRS is compared.
The second case is reproduction of the result of \cite{lee97b}, who investigated Raman scattering
of O VI in symbiotic stars.
In this case, the H I region is an expanding spherical wind around the giant component.
Schematic illustrations of the two cases are shown in  Fig.~\ref{fig:shell_model}.

\subsection{Formation of Balmer Wings in a Static Spherical H~I Region}

The central point-like far UV source surrounded by a spherical H I region with radius $R$ is characterized by a flat continuum.
The static neutral region is assumed to be of
uniform H I density $n_{\rm HI}$. We fix the radial column density $N_{\rm HI}$ defined by
\begin{equation}
N_{\rm HI}=n_{\rm HI} R = 10^{23}{\rm\ cm^{-2}},
\end{equation}
and vary the number of cells. We set $N_{xyz}=N_x=N_y=N_z$,
so that the total number of cells is given by $N_{xyz}^3$.
We generate $10^7$ photons for each simulation.
The initial photons are generated at the center of the H I sphere in accordance with
	the emisivity $j_{e}$ given by the three dimensional Dirac delta function
	\begin{equation}
j_{e}(\lambda,x,y,z) =
{\lambda \over {\lambda_{max}-\lambda_{min}}} \delta^{(3)}({\bf r}),
\end{equation}
where $\lambda_{max}$ and $\lambda_{min}$ are the maximum and minmum wavelengths of initial photons.

Fig.~\ref{fig:shell_spec} shows optical spectra formed
through Raman scattering for the cases of $N_{xyz}=3, 10$ and 100.
The vertical axis shows the Raman conversion efficiency ($RCE$)
defined as the number ratio per unit wavelength of incident far UV photons and Raman-scattered optical photons.
The left and right panels show Balmer wings formed through Raman scattering
around H$\alpha$ and H$\beta$, respectively.
The solid lines represent the analytic solutions and open circles
show simulation results obtained using STaRS.
The analytic solutions are obtained from the non-grid based simulation in \cite{chang15,chang18}.

Far from the line centers, the simulation results for $N_{xyz}=3$ are slightly higher
than the analytic solutions.
Other than this, the agreement is fairly good,
indicating little dependence on $N_{xyz}$ in the case of a static H I region.
The volume of the H I region with $N_{xyz}=3$ is larger than the sphere with the radius $R$
as the H~I density associated with a cell is determined by the central position $(X_C,Y_C,Z_C)$ of the cell.
When $(X_C^2 +Y_C^2 +Z_C^2)^{1/2}$ is smaller than $R$, the H~I density of the cell is assigned to be $n_{HI}$.
Otherwise, the density is assigned to be zero.
The number ratio between Raman scattered photons near H$\alpha$ and the total initial photons near Ly$\beta$ is
20.77 \% for the analytic solution, whereas the simulations give 21.17, 20.91, and  20.77 \% for $N_{xyz} = 3$, 10, and 100,
respectively.

In Fig.~\ref{fig:image_ha}, we present the polarimetric data of Raman scattered H$\alpha$ projected to the celestial sphere.
The surface brightness, the degree of polarization $p$, and the direction of polarization are shown in the left, middle and the right panels, respectively.
Here, in terms of the Stokes parameters,
the degree of polarization $p$ and the position angle $\phi_p$ are given by
\begin{equation}
    p = \frac{\sqrt{Q^2 +U^2}}{I}, \quad \phi_p = \frac{1}{2} \tan^{-1}\left( \frac{U}{Q} \right).
\end{equation}
In the left panels, the surface brightness is shown in logarithmic scale.
With our choice of a rather large value of $N_{\rm HI}$,
a considerable fraction of Raman-scattered photons are formed near the source, which
leads to excellent agreement between the analytic results and those obtained using STaRS.
One may notice increase in smoothness of the surface brightness as $N_{xyz}$ increases.
It is also noticeable that $p$ becomes large with increasing distance from the center.
The concentric polarization pattern reflects the spherically symmetric scattering geometry.
In this particular case, $N_{xyz}>10$ appears to be sufficient to describe the the analytic result of the static medium.

\begin{figure*}
	\includegraphics[width=170mm]{Model_Scheme.pdf}
	\caption{Schematic illustrations of two cases to test the code.
	}
	\label{fig:shell_model}
\end{figure*}

\begin{figure*}
	\includegraphics[width=170mm]{spec.pdf}
	\caption{Raman conversion efficiency of H$\alpha$ (left panel) and H$\beta$ (right panel) broad wings of the static spherical case.
	The black solid lines are the spectra by the analytic method.
	The $x$ axis is the wavelength.
	The red, blue, and green open circles are the spectra by STaRS for $N_{xyz}= 3$, 10, and 100.
	}
	\label{fig:shell_spec}
\end{figure*}

\subsection{Raman Scattering of O VI in Expanding H I Region}

\cite{lee97b} presented their basic study of line formation of Raman O VI  in a symbiotic star consisting of a white dwarf and a mass losing giant.
The O VI emission region
near the white dwarf component was assumed to be a point source broadened thermally with $T=10^4\, \rm K$.
For simplicity, the slow stellar wind from the red giant component is assumed to be entirely neutral ignoring
the photoionization by the white dwarf.
The orbital separation is $10R_*$ of which $R_*$ is the radius of the red giant.
The positions of the white dwarf and the red giant are $(0,0,10R_*)$ and $(0,0,0)$, respectively.
The initial photons are generated at the position of the white dwarf according to
	the emissivity $j_{e}$ given by
	\begin{eqnarray}
		&&j_{e}(\lambda,x,y,z) = \\
		\nonumber
		&&{1 \over {\sigma_{th} \sqrt{2\pi}}}  \exp\left[{(\lambda - \lambda_{1032})^2 \over {2 \sigma_{th}^2}} \right] \delta^{(3)}(x,y,z-10R_*)
	\end{eqnarray}
	where $\lambda_{1032}$ is the center wavelength of O VI $\lambda$ 1032 and
	$\sigma_{th} \sim 0.008{\rm\ \AA}$ is the thermal width of O VI $\lambda$ 1032 with $T = 10^4 K$.

The velocity $\bf v(r)$ and H I number density $n({\bf r})$ are the functions of the distance from the red giant $r=|\bf r|$
 are given by
\begin{eqnarray}
    {\bf v}({\bf r}) &=& v_{\infty}(1 - R_*/r) \frac{\bf r}{ r}\\
    \nonumber
    \quad
    n({\bf r}) &=& n_0 (R_*/r)^{2} (1 - R_*/r)^{-1},
\end{eqnarray}
where $n_0$ is the characteristic number density defined in Eq.~2.11 of \cite{lee97b}
and $v_{\infty}$ is the terminal velocity.
As a second code test case, we revisit Raman O VI formation illustrated in Figs.~1 and 6 of \cite{lee97b}.
Fixing $v_\infty=20\, \rm km\, s^{-1}$, we consider the three values of
$\tau_0 = 0.5 $, 1, and 10, where $\tau_0 = n_0 R_* \sigma_{\rm tot}$.
We generate $10^9$ photons for each simulation.

Fig.~\ref{fig:lee_1997} shows the spectra and the Stokes parameter
$P=Q/I$ for Raman scattered O VI at 6825 \AA.
It should be noted that the integrated $U$ vanishes due to the axial symmetry about the $z$ axis.
Therefore, the signed ratio $Q/I$ represents the degree and direction of polarization, where
a positive  and a negative $Q$ correspond to  the polarization in direction perpendicular and parallel to the $z$ axis, respectively.
The observer's line of sight lies in $x$-$y$ plane.
Noting that it is perpendicular to the symmetry axis,
maxium polarization can be developed along this direction or the symmetry $z-$axis.
We collecte the photons escaping toward the observer.
The fraction of the collected Raman photons near 6825 \AA\, is
0.11, 0.21, and 1.82 \% for $\tau_0 = 0.5$, 1, and 10, respectively.

The line profiles obtained using STaRS differ slightly from those \cite{lee97b} presented.
In particular,
the red peaks obtained from STaRS are more enhanced than those presented by \cite{lee97b},
which is attributed to the geometrical truncation adopted by them.
The difference is in the range of 10-15 per cent with respect to the red peak.
We find that the agreement gets better where the full range of the scattering region is taken into account.
In the bottom panels, we find overall agreement in the polarization behaviors.
The noisy features with $\Delta V > 200 \rm\, km\, s^{-1}$ are attributed to the small number statistics of
collected photons.
Blue photons are scattered mostly in the compact and dense region between the red giant and white dwarf, resulting in
development of strong polarization in the direction perpendicular to the $z$ axis.
In contrast, red photons are scattered mainly in a quite extended region near the red giant,
leading to weak polarization and enhanced line flux.

Fair agreement shown in  Fig.~\ref{fig:shell_spec} and \ref{fig:lee_1997} demonstrates
that STaRS has been well-tested. Furthermore, as illustrated in
Fig.~\ref{fig:image_ha}, STaRS is capable of study of radiative transfer for spectropolarimetric imaging observations.

\begin{figure*}
	\includegraphics[width=165mm]{image_Ha.pdf}
	\caption{The surface brightness (left), the degree of polarization (center), and the direction of polarization (right) of the projected H$\alpha$ photons.
    The panels in the first, second, and third rows represent the results obtained using STaRS for $N_{xyz}= 3$, 10, and 100. The bottom panels represent the result obtained using an analytic method.
	}
	\label{fig:image_ha}
\end{figure*}

\begin{figure*}
	\includegraphics[width=170mm]{com_over_6825.eps}
	\caption{Line formation of Raman O VI 6825 features in an expanding H I region.
	The blue solid lines are computed by 'STaRS'.
	The red dashed lines are the results of Fig.~6 in \cite{lee97b}.
	}
	\label{fig:lee_1997}
\end{figure*}

\section{Summary}

We have developed a 3D grid-based Monte Carlo code 'STaRS' for radiative transfer through Raman
and Rayleigh scattering, which can be mainly used to investigate line formation of Raman-scattered
features in a thick neutral region illuminated by a strong far UV emission source. Favorable
conditions for Raman scattering with atomic hydrogen are easily met in symbiotic stars, young planetary nebulae and active galactic nuclei.
Through a couple of tests, we have successfully demonstrated that 'STaRS' is a flexible
code to deal with radiative transfer in a thick neutral media yielding multidimensional
spectropolarimetric and imaging data.
'STaRS' is easily accessed in Github 'https://github.com/csj607/STaRS'.

\acknowledgments

This research was supported by the Korea Astronomy and Space Science Institute
under the R\&D program (Project No. 2018-1-860-00) super-vised
by the Ministry of Science, ICT and Future Planning.
This work was also supported by a National Research Foundation of Korea (NRF) grant
funded by the Korea government (MSIT; No. NRF-2018R1D1A1B07043944).
Seok-Jun is very grateful to Dr. Kwang-Il Seon for his help
in adoption of grid-based technique.