# End-to-End Multi-Speaker Speech Recognition using Speaker Embeddings and
Transfer Learning

###### Abstract

This paper presents our latest investigation on end-to-end automatic speech
recognition (ASR) for overlapped speech. We propose to train an end-to-end
system conditioned on speaker embeddings and further improved by transfer
learning from clean speech. This proposed framework does not require any
parallel non-overlapped speech materials and is independent of the number of
speakers. Our experimental results on overlapped speech datasets show that
joint conditioning on speaker embeddings and transfer learning significantly
improves the ASR performance.

Index Terms: end-to-end asr, overlapped speech

##  1 Introduction

Recently, deep learning technology has boosted automatic speech recognition
(ASR) performance significantly [1, 2, 3, 4]. Overlapped speech – well known
in a more general context as the cocktail party problem – remains, however, to
be a largely unsolved problem. Its difficulty can be mainly explained by high
similarity of acoustic characteristics of signals that need to be separated
and absence of any other obvious clues that might guide a potential solution.

Speech recognition of overlapped speech is usually approached in two stages:
first overlapped speech is split to separate recordings of each speaker, then
speech recognition is performed on separated recordings. Previous works on
speech separation problem include computational auditory scene analysis [5],
non-negative matrix factorization [6], graphical modeling [7] and spectral
clustering [8]. Similarly to speech recognition, speech separation methods
have also made major progress with the help of deep learning. The deep
clustering method has been introduced in [9] and consequently improved in [10,
11] and has become one of the most remarkable speech separation methods in the
recent years. Deep clustering operates on a spectrogram of overlapped speech
and learns to map time-frequency (T-F) units to high dimensional embedding
space. More recently, a different approach, named VoiceFilter, has been
proposed in [12]. The work simplifies the problem of multiclass classification
of T-F units over multiple speakers to binary classification between target
speaker’s speech and everything else. To condition the neural network on
specific speaker, the input is extended with a speaker embedding vector, which
is extracted by a separately trained network from the reference clean speech.

A new type of ASR systems, called end-to-end ASR, has emerged in the past
years [13, 14, 15, 16, 17]. End-to-end ASR maps the acoustic signal to written
language with a single encoder-decoder recurrent neural network and does not
require any domain specific knowledge for solving of intermediate subtasks,
such as grapheme-to-phoneme conversion. Recently, two works have proposed to
integrate speaker separation stage to end-to-end ASR. The first one [18]
connects pretrained deep clustering model and end-to-end ASR for the
subsequent join fine-tuning for the better ASR results. The second one [19]
removes explicit speech separation part and trains end-to-end ASR for
simultaneous speech separation and recognition by permutation invariant
procedure, in which ASR system is optimized for multiple outputs corresponding
to multiple speakers in the input mixture. Joint speech separation and
recognition is also described in [20], although this work is based on the
conventional ASR. Both [19] and [20] suggest that transfer learning from clean
speech ASR improves the results of overlapped speech ASR.

Our work blends the ideas from [21, 22, 12, 18, 19, 20] and proposes to train
an end-to-end overlapped speech recognition system conditioned by speaker
embeddings and improved by transfer learning from clean speech. It has
advantages over [12] [18] and [20] of not requiring parallel clean speech
material and over [19] of not depending on the number of speakers. We evaluate
our proposed framework on overlapped speech datasets with two and three
overlapped speakers, within and across settings. Overall, we observe
significant improvements over the baseline end-to-end ASR system.

##  2 Method

The outline of the proposed method is presented on Figure 1. It shows two
separate neural network models, the speaker encoder and the end-to-end ASR,
together with their inputs and outputs. The speaker encoder takes reference
speech utterances on the input and produces speaker embedding vectors, as
described in sections 2.2 and 3.3. The end-to-end ASR takes acoustic features
of overlapped speech with speaker embedding vector of the target speaker on
the input and generates transcription of the target speaker’s speech on the
output, which is used to update the parameters of end-to-end ASR model during
the training or provided as the final output during the decoding. A more
detailed description of end-to-end ASR is given in sections 2.1 and 3.2.
Optionally, clean speech can be used during the training in order to perform a
basic form of transfer learning, which is described in section 2.3.

Figure 1:  Overview of end-to-end ASR using speaker embeddings and transfer
learning.

###  2.1 End-to-end ASR

End-to-end ASR has a hybrid CTC/attention architecture, described in detail in
[17]. The input of the model is defined as a $T$-length sequence of $d$
dimensional feature vectors $X=\\{x_{t}\in\mathbb{R}^{d}|t=1,\dots,T\\}$, and
the output of the model is defined as a $U$-length sequence of output labels
$Y=\\{y_{u}\in\mathcal{U}|u=1,\dots,U\\}$, where $\mathcal{U}$ is a set of
distinct output labels and usually $U<T$. During the training, a weighted sum
of CTC loss and attention-based cross-entropy loss is minimized:

$\displaystyle\mathcal{L}=\lambda\mathcal{L}_{\text{ctc}}+(1-\lambda)\mathcal{L}_{\text{att}},$
(1)

where $0\leq\lambda\leq 1$. Attention-based cross-entropy loss is calculated
according to predictions of the attention-based encoder-decoder network:

$\displaystyle\mathcal{L}_{\text{att}}=$
$\displaystyle-\log\mathnormal{p}_{\text{att}}(Y|X)$ (2)
$\displaystyle\mathnormal{p}_{\text{att}}(Y|X)=$
$\displaystyle\prod_{u}\mathnormal{p}(y_{u}|X,y_{1:u-1})$ (3)
$\displaystyle\mathnormal{p}(y_{u}|X,y_{1:u-1})=$ $\displaystyle\text{
Decoder}(\mathbf{r}_{u},\mathbf{q}_{u-1},y_{u-1})$ (4)
$\displaystyle\mathbf{h}_{t}=$ $\displaystyle\text{ Encoder}(X)$ (5)
$\displaystyle a_{ut}=$ $\displaystyle\text{
Attention}(\\{a_{u-1}\\}_{t},\mathbf{q}_{u-1},\mathbf{h}_{t})$ (6)
$\displaystyle\mathbf{r}_{u}=$ $\displaystyle\sum_{t}a_{ut}\mathbf{h}_{t}.$
(7)

Here, $\text{Encoder}(\cdot)$ and $\text{Decoder}(\cdot)$ are recurrent neural
networks, $\text{Attention}(\cdot)$ is an attention mechanism and
$\mathbf{h}_{t}$, $\mathbf{q}_{u-1}$ and $\mathbf{r}_{u}$ are the hidden
vectors. CTC loss is calculated from a linear transformation of encoder output
and all possible $T$-length sequences of an extended output labels set
$Z=\\{z_{t}\in\mathcal{U}\cup\texttt{<blank>}|t=1,\dots,T\\}$, corresponding
to sequence $Y$ of original output labels:

$\displaystyle\mathcal{L}_{\text{ctc}}$
$\displaystyle=-\log\mathnormal{p}_{\text{ctc}}(Y|X)$ (8)
$\displaystyle\mathnormal{p}(Y|X)$
$\displaystyle\approx\underbrace{\sum_{Z}\prod_{t}\mathnormal{p}(z_{t}|z_{t-1},Y)\mathnormal{p}(z_{t}|X)}_{\triangleq\mathnormal{p}_{\text{ctc}}(Y|X)}\mathnormal{p}(Y)$
(9) $\displaystyle\mathnormal{p}(z_{t}|X)$ $\displaystyle=\text{
Softmax}(\text{Lin}(\mathbf{h}_{t})).$ (10)

During the decoding, the same CTC and attention-based probabilities are summed
with possibly another weight and are used to find the most probable output
labels sequence:

$\displaystyle\hat{Y}=\operatorname*{arg\,max}_{Y}\\{\lambda\log\mathnormal{p}_{\text{ctc}}(Y|X)+(1-\lambda)\log\mathnormal{p}_{\text{att}}(Y|X)\\}.$
(11)

###  2.2 Speaker embeddings

Speaker embedding is a vector of fixed dimensionality that represents
speaker’s characteristics and can be extracted from a reference recording of
speaker’s speech. Speaker embeddings have been shown to be a useful source of
information about speaker in many tasks, including speaker verification,
speaker diarization [23], speech synthesis [24] and speech separation [12]. We
condition ASR system for the recognition of speech of a certain speaker in the
recording of overlapped speech. This arrangement removes a major roadblock of
the permutation problem, appearing when multiple correct outputs are possible
for a single input, and simplifies application of a wide range of well studied
machine learning methods for further system optimization.

###  2.3 Transfer learning

Speech recognition of overlapped speech can be viewed as speech recognition of
clean speech in mismatched conditions and be addressed by numerous transfer
learning methods. It has been proposed in [20] to utilize teacher-student
training for transfer learning from clean to overlapped speech recognition.
This approach has a limitation of applicability on training sets with parallel
clean and overlapped speech only, which is hardly achievable in real life
scenarios. Another example of transfer learning from clean to overlapped
speech recognition is given in [19], which applies parameters transfer. Multi-
condition training is an alternative to the parameters transfer method, in
which training samples from different conditions are mixed together and used
for the training simultaneously. Multi-condition training has been
demonstrated to provide better results for transfer learning between languages
[25] as well as between channels [26]. In this work, we experiment with both
parameters transfer and multi-condition training methods in order to improve
overlapped speech recognition by utilization of clean speech training data.

##  3 Experimental setup

###  3.1 Datasets

We evaluate our models on the widely used mixed speech datasets wsj0-2mix and
wsj0-3mix [9, 10]. The datasets contain two-speaker and three-speaker mixtures
of randomly selected utterances from WSJ0 corpus. Training, development and
evaluation sets, named tr, cv and tt, are generated from WSJ0 training,
development and evaluation sets si_tr_s, si_dt_05 and si_et_05 and contain
speech of speakers of both genders in different combinations. We use max
version of the datasets, meaning that the length of every mixed speech
utterance is chosen to be maximum of the lengths of original utterances being
used for a mixture. The sampling rate of the used datasets is 16 kHz.
Training, development and evaluation sets of both datasets contain 20000, 5000
and 3000 utterances. Training and development sets contain speech of same 101
speakers, therefore development sets are used for the evaluation in closed
speaker set condition. Evaluation sets contain speech of another 19 speakers,
therefore evaluation sets are used for the evaluation in open speaker set
condition. Total durations of training, development and evaluation sets of
wsj0-2mix dataset are 46, 11 and 7 hours. Total durations of training,
development and evaluation sets of wsj0-3mix dataset are 51, 13 and 8 hours.
LibriSpeech [27] train-clean-100 dataset is used for the transfer learning
experiments. It contains 28539 utterances of read speech by 251 speakers and
has total duration of 100 hours. The sampling rate of the dataset is 16 kHz.

###  3.2 Baseline

Our end-to-end ASR system is based on ESPnet toolkit [28] and its WSJ recipe.
80-dimensional log Mel filterbank coefficients with pitch with a frame length
of 25 ms and shift of 10 ms are used as input features. The input features are
extracted and normalized to zero mean and unit variance with Kaldi toolkit
[29]. The encoder network consists of four BLSTM layers with 320 units in each
layer and direction and linear projection layer with 320 units. No subsampling
is applied to the input. The decoder network consists of one LSTM layer with
300 units. Additive attention mechanism with 320 dimensions is utilized. We
use 49 characters as output units. PyTorch backend of ESPnet is used to
implement the networks. Training is performed with AdaDelta optimizer [30] and
gradient clipping on two GPUs in parallel with a batch size of 30 for 30
epochs. The optimizer is initialized with $\rho=0.95$ and $\epsilon=10^{-8}$.
$\epsilon$ is halved after an epoch if performance of the model did not
improve on development set. The model with the highest accuracy on development
set is used for the decoding. The CTC weight $\lambda$ is set to $0.2$ during
the training and $0.3$ during the decoding. The decoding is performed with a
beam search with a beam size of 30.

The decoding makes use of word-based RNN-LM [31] with a weight of $1.0$. Word-
based RNN-LM is trained on the LM training set of WSJ0 corpus containing 37M
words and 1.6M sentences, dictionary size is 65K words. Word-based RNN-LM
contains 1 LSTM layer with 1000 units. The stochastic gradient descent
optimizer is used to train word-based RNN-LM with a batch size of 300 for 20
epochs.

###  3.3 Speaker embeddings extractor

We extract 512-dimensional speaker embeddings from the reference utterances
with x-vector system from Kaldi toolkit. We use the pretrained model
downloaded from http://kaldi-asr.org/models/m8. The model is trained on
augmented VoxCeleb 1 [32] and VoxCeleb 2 [33] datasets according to the
procedure closely following the description in [34] and evaluated on Speakers
in the Wild dataset [35] with 3.5% equal error rate. The input features of
x-vector extractor are 30-dimensional MFCCs without cepstral truncation with a
frame length of 25 ms and shift of 10 ms. Mean normalization with a sliding
window of up to three seconds is applied to the input features. Speaker
embeddings are extracted from voiced frames only, which are selected by the
same energy-based VAD. We obtain one vector from each utterance and average
them per speaker to get speaker specific vector. $L_{2}$ normalized speaker
embedding is used as an additional input for the ASR model. We denote
insertion of the speaker embedding vector to the beginning of the sequence of
acoustic features vectors as horizontal stacking. We denote concatenation of
the speaker embedding vector with every acoustic features vector as vertical
stacking. Horizontal stacking requires speaker embedding and acoustic features
to have same size, which can be achieved either by downscaling of the speaker
embedding vector, which we perform with the trainable linear transformation,
or by padding of the acoustic feature vectors, which we perform by appending
appropriate number of zeros to each acoustic features vector. Vertical
stacking allows the size of the speaker embedding to be independent of the
size of the acoustic features.

##  4 Results

###  4.1 Speaker embeddings inclusion strategies

The first set of experiments aims to determine the best strategy for inclusion
of speaker embeddings in the model’s input. While vertical stacking does not
enforce same size of the speaker embedding and acoustic features, we perform
two experiments with vertical stacking: the first one with unchanged sizes of
the input vectors and the second one with the downscaled speaker embedding.
The second experiment isolates the effect of different stacking types from the
effect of the speaker embedding downscaling. 50 utterances of reference speech
are used to produce each speaker embedding in this experiment. It can be seen
from the data in Table 1 that vertical stacking clearly outperforms horizontal
stacking, while the downscaling of the speaker embedding results in minor
degradation of the model’s performance.

Table 1:  Results (WER, %) with different speaker embedding inclusion
strategies on the two-speaker overlapped speech dataset.
---
Strategy | dev | eval
Baseline (no speaker embeddings) | 79.6 | 85.7
Horizontal stacking with downscaled embedding | 77.0 | 84.1
Horizontal stacking with padded acoustic features | 83.0 | 88.6
Vertical stacking with downscaled embedding | 11.7 | 24.9
Vertical stacking with unchanged vectors’ sizes | 11.4 | 22.1

###  4.2 Amount of reference speech data

The next set of experiments is concerned with amount of reference speech data
required for the generation of speaker embeddings. Table 2 gives an overview
of ASR performance with different numbers of reference utterances in case of
two-speaker overlapped speech. It is apparent from this table that larger
amount of reference speech data allows to generate more general speaker
embeddings preventing the ASR model from the overfitting towards known
speakers and seen utterances. It is worth noting, however, that even one
reference utterance of the duration of approximately ten seconds appears to be
sufficient material for the major improvement.

Table 2:  Results (WER, %) with different amount of reference speech data on
the two-speaker overlapped speech dataset.
---
Reference speech amount per speaker | dev | eval
Utterances | Seconds | Voiced frames |  |
1 | 8.3 ±2.9 | 641 ±271 | 15.8 | 32.6
5 | 42.3 ±6.4 | 3156 ±623 | 17.3 | 29.1
10 | 84.6 ±10.1 | 6341 ±1090 | 11.3 | 22.6
20 | 170.3 ±16.9 | 12545 ±1783 | 10.8 | 22.5
50 | 426.3 ±34.2 | 31490 ±4260 | 11.4 | 22.1

###  4.3 Transfer learning

We evaluate two transfer learning approaches, namely parameters transfer and
multi-condition training, for improving ASR performance on overlapped speech
by utilization of non-parallel clean speech training data. Table 3 presents
the results of the systems trained on training sets of wsj0-2mix and wsj-3mix
datasets and tested on development and evaluation sets of the corresponding
datasets. The most striking observation to emerge from the results is that the
training process has not converged on the dataset with three overlapping
speakers, but the system has been able to decode the same recordings when
trained on the combination of overlapped and clean speech datasets. This
finding can be attributed to higher complexity of the modeled function in case
of increased number of overlapping speakers, which the neural network could
not learn just from the overlapped speech data, and demonstrates how crucial
the role of the transfer learning approach for the solution of certain
problems can be. The transfer learning results are expectedly better on the
dataset with two overlapping speakers as well, especially for the open speaker
condition, what is due to the number of additional speaker embeddings in the
training data and subsequent better generalization of the relationship between
speaker embedding and relevant acoustic features. Finally, multi-condition
training has demonstrated slightly better results than parameters transfer.
This finding is in agreement with the previous reports on transfer learning
applications in ASR.

Table 3:  Results (WER, %) of the baseline ASR, speaker embeddings
conditioning and transfer learning on the two- and three-speaker overlapped
speech datasets.
---
| wsj0-2mix | wsj0-3mix
| dev | eval | dev | eval
Baseline | 79.6 | 85.7 | 95.9 | 96.0
\+ speaker embeddings | 11.4 | 22.1 | 95.6 | 95.7
\+ parameters transfer | 8.8 | 16.9 | 22.7 | 45.3
\+ multi-condition training | 8.5 | 14.6 | 21.7 | 42.9

As our method does not utilize any explicit knowledge about the number of
overlapping speakers, the models should theoretically also work for testing
data with larger or smaller number of overlapping speakers than in the
training data. We test whether this is true in practice by decoding wsj0-3mix
testing data with the model trained on wsj0-2mix training data (combined with
LibriSpeech 100) and vice versa. Encouraged by the success of the previous
transfer learning experiments, we also train a system on a combination of
clean and overlapped two- and three-speaker datasets. The results of testing
with mismatching number of overlapping speakers are given in Table 4. In
general it seems that the proposed method does not depend on the number of
overlapping speakers but can benefit from training on larger amount of
different conditions of speech overlap. A possible explanation of the slightly
worse result of the best system on the open speaker set condition with two
overlapping speakers might be a bias towards WSJ0 speakers in the combined
training dataset introduced by the addition of wsj0-3mix dataset, which is
slightly larger than wsj0-2mix dataset and therefore has more impact on two-
speaker testing data than wsj0-2mix on three-speaker testing data in this
experiment.

Table 4:  Results (WER, %) of the proposed method depending on the training
data.
---
Training data | wsj0-2mix | wsj0-3mix
| dev | eval | dev | eval
LibriSpeech 100 + wsj0-2mix | 8.5 | 14.6 | 45.2 | 55.3
LibriSpeech 100 + wsj0-3mix | 7.8 | 28.5 | 21.7 | 42.9
LibriSpeech 100 + wsj0-{2,3}mix | 4.8 | 15.2 | 15.5 | 32.3

###  4.4 Comparison with earlier work

Table 5 compares our best result on the evaluation set of wsj0-2mix dataset
with the results reported on the same dataset in the previous works. From the
table we can see that the proposed system outperforms the best known result by
42% relatively. However, it should be noted that the listed systems differ
from each other in a number of ways, including types of ASR system, types of
LM and types and amount of training data.

Table 5:  Results (WER, %) of the proposed method and previous works on the
two-speaker overlapped speech dataset.
---
Method | eval
Deep Clustering, hybrid ASR [10] | 30.8
Permutation Invariant Training, hybrid ASR [36] | 28.2
Permutation Invariant Training, end-to-end ASR [19] | 28.2
Speaker Parallel Attention, end-to-end ASR [37] | 25.4
Proposed, end-to-end ASR | 14.6

(a) Without speaker embedding

(b) With the speaker embeddings for the speaker 01z

(c) With the speaker embeddings for the speaker 20h

Figure 2:  Visualization of the hidden vector sequences for the utterance
01zc020o_2.3474_20hc010j_-2.3474 of wsj-2mix dataset.

Following [19], we also visualize the encoder networks outputs of an example
utterance (see Figure 2). We apply principal component analysis to the hidden
vectors on the vertical axis as well. Figure 2(a) shows the output of the
baseline’s encoder, while figures 2(b) and 2(c) show the encoder network
conditioned on speaker embeddings of two different speakers. Some patterns
from the encoder’s output of the baseline model appear on the conditioned
encoder’s output for the first speaker, and another ones appear on the
encoder’s output for the second speaker. This observation suggests that the
conditioning on the speaker embeddings indeed allows the encoder network to
perform the separation of overlapped speech.

##  5 Conclusions

In this paper, we proposed an effective end-to-end speech recognition
framework for overlapped speech using speaker embeddings and transfer learning
techniques. Experimental results on simulated overlapped speech datasets
revealed that using speaker embeddings our framework was able to automatically
identify relevant information of the target speaker for recognition. The
application of transfer learning technique played a crucial role with the
increasing number of speakers. Finally, we observed significant improvements
over the baseline end-to-end system even while using just ten seconds of
reference speech per speaker.

## References

  * [1] G. Hinton, L. Deng, D. Yu, G. Dahl, A. rahman Mohamed, N. Jaitly, A. Senior, V. Vanhoucke, P. Nguyen, T. Sainath, and B. Kingsbury, “Deep Neural Networks for Acoustic Modeling in Speech Recognition,” _Signal Processing Magazine_ , 2012.
  * [2] G. Dahl, D. Yu, L. Deng, and A. Acero, “Context-Dependent Pre-trained Deep Neural Networks for Large Vocabulary Speech Recognition,” in _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 20, 2012, pp. 30–42.
  * [3] W. Xiong, J. Droppo, X. Huang, F. Seide, M. Seltzer, A. Stolcke, D. Yu, and G. Zweig, “The Microsoft 2016 Conversational Speech Recognition System,” _CoRR_ , 2016.
  * [4] D. Povey, V. Peddinti, D. Galvez, P. Ghahremani, V. Manohar, X. Na, Y. Wang, and S. Khudanpur, “Purely Sequence-Trained Neural Networks for ASR Based on Lattice-Free MMI,” in _Interspeech_ , 2016, pp. 2751–2755.
  * [5] D. Wang and G. J. Brown, _Computational auditory scene analysis: Principles, algorithms, and applications_. Wiley-IEEE press, 2006.
  * [6] M. N. Schmidt and R. K. Olsson, “Single-channel speech separation using sparse non-negative matrix factorization,” in _Ninth International Conference on Spoken Language Processing_ , 2006.
  * [7] J. R. Hershey, S. J. Rennie, P. A. Olsen, and T. T. Kristjansson, “Super-human multi-talker speech recognition: A graphical modeling approach,” _Computer Speech & Language_, vol. 24, no. 1, pp. 45–66, 2010.
  * [8] F. R. Bach and M. I. Jordan, “Learning spectral clustering, with application to speech separation,” _Journal of Machine Learning Research_ , vol. 7, no. Oct, pp. 1963–2001, 2006.
  * [9] J. R. Hershey, Z. Chen, J. Le Roux, and S. Watanabe, “Deep clustering: Discriminative embeddings for segmentation and separation,” in _Acoustics, Speech and Signal Processing (ICASSP), 2016 IEEE International Conference on_. IEEE, 2016, pp. 31–35.
  * [10] Y. Isik, J. Le Roux, Z. Chen, S. Watanabe, and J. R. Hershey, “Single-Channel Multi-Speaker Separation Using Deep Clustering,” _Interspeech 2016_ , pp. 545–549, 2016.
  * [11] Z.-Q. Wang, J. Le Roux, and J. R. Hershey, “Alternative Objective Functions for Deep Clustering,” in _Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2018.
  * [12] Q. Wang, H. Muckenhirn, K. Wilson, P. Sridhar, Z. Wu, J. Hershey, R. A. Saurous, R. J. Weiss, Y. Jia, and I. L. Moreno, “VoiceFilter: Targeted Voice Separation by Speaker-Conditioned Spectrogram Masking,” _arXiv preprint arXiv:1810.04826_ , 2018.
  * [13] A. Hannun, C. Case, J. Casper, B. Catanzaro, G. Diamos, E. Elsen, R. Prenger, S. Satheesh, S. Sengupta, A. Coates _et al._ , “Deep speech: Scaling up end-to-end speech recognition,” _arXiv preprint arXiv:1412.5567_ , 2014\.
  * [14] Y. Miao, M. Gowayyed, and F. Metze, “EESEN: End-to-end speech recognition using deep RNN models and WFST-based decoding,” in _Automatic Speech Recognition and Understanding (ASRU), 2015 IEEE Workshop on_. IEEE, 2015, pp. 167–174.
  * [15] D. Bahdanau, J. Chorowski, D. Serdyuk, P. Brakel, and Y. Bengio, “End-to-end attention-based large vocabulary speech recognition,” in _Acoustics, Speech and Signal Processing (ICASSP), 2016 IEEE International Conference on_. IEEE, 2016, pp. 4945–4949.
  * [16] W. Chan, N. Jaitly, Q. Le, and O. Vinyals, “Listen, attend and spell: A neural network for large vocabulary conversational speech recognition,” in _Acoustics, Speech and Signal Processing (ICASSP), 2016 IEEE International Conference on_. IEEE, 2016, pp. 4960–4964.
  * [17] S. Watanabe, T. Hori, S. Kim, J. R. Hershey, and T. Hayashi, “Hybrid CTC/attention architecture for end-to-end speech recognition,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 11, no. 8, pp. 1240–1253, 2017.
  * [18] S. Settle, J. Le Roux, T. Hori, S. Watanabe, and J. R. Hershey, “End-to-end multi-speaker speech recognition,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2018, pp. 4819–4823.
  * [19] H. Seki, T. Hori, S. Watanabe, J. Le Roux, and J. R. Hershey, “A Purely End-to-End System for Multi-speaker Speech Recognition,” in _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , vol. 1, 2018, pp. 2620–2630.
  * [20] Z. Chen, J. Droppo, J. Li, and W. Xiong, “Progressive joint modeling in unsupervised single-channel overlapped speech recognition,” _IEEE/ACM Transactions on Audio, Speech and Language Processing (TASLP)_ , vol. 26, no. 1, pp. 184–196, 2018.
  * [21] K. Žmolíková, M. Delcroix, K. Kinoshita, T. Higuchi, A. Ogawa, and T. Nakatani, “Speaker-Aware Neural Network Based Beamformer for Speaker Extraction in Speech Mixtures,” 2017, pp. 2655–2659.
  * [22] J. Wang, J. Chen, D. Su, L. Chen, M. Yu, Y. Qian, and D. Yu, “Deep Extractor Network for Target Speaker Recovery from Single Channel Speech Mixtures,” _Proc. Interspeech 2018_ , pp. 307–311, 2018.
  * [23] Q. Wang, C. Downey, L. Wan, P. A. Mansfield, and I. L. Moreno, “Speaker diarization with LSTM,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 5239–5243.
  * [24] Y. Jia, Y. Zhang, R. Weiss, Q. Wang, J. Shen, F. Ren, P. Nguyen, R. Pang, I. L. Moreno, Y. Wu _et al._ , “Transfer learning from speaker verification to multispeaker text-to-speech synthesis,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 4485–4495.
  * [25] G. Heigold, V. Vanhoucke, A. Senior, P. Nguyen, M. Ranzato, M. Devin, and J. Dean, “Multilingual acoustic models using distributed deep neural networks,” in _ICASSP_. IEEE, 2013, pp. 8619–8623.
  * [26] P. Ghahremani, V. Manohar, H. Hadian, D. Povey, and S. Khudanpur, “Investigation of transfer learning for ASR using LF-MMI trained neural networks,” in _Automatic Speech Recognition and Understanding Workshop (ASRU), 2017 IEEE_. IEEE, 2017, pp. 279–286.
  * [27] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, “Librispeech: an ASR corpus based on public domain audio books,” in _Acoustics, Speech and Signal Processing (ICASSP), 2015 IEEE International Conference on_. IEEE, 2015, pp. 5206–5210.
  * [28] S. Watanabe, T. Hori, S. Karita, T. Hayashi, J. Nishitoba, Y. Unno, N.-E. Y. Soplin, J. Heymann, M. Wiesner, N. Chen _et al._ , “ESPnet: End-to-End Speech Processing Toolkit,” _Proc. Interspeech 2018_ , pp. 2207–2211, 2018.
  * [29] D. Povey, A. Ghoshal, G. Boulianne, L. Burget, O. Glembek, N. Goel, M. Hannemann, P. Motlicek, Y. Qian, P. Schwarz _et al._ , “The Kaldi speech recognition toolkit,” in _ASRU_ , no. EPFL-CONF-192584. IEEE Signal Processing Society, 2011.
  * [30] M. D. Zeiler, “Adadelta: an adaptive learning rate method,” _arXiv preprint arXiv:1212.5701_ , 2012.
  * [31] T. Hori, J. Cho, and S. Watanabe, “End-to-end Speech Recognition with Word-based RNN Language Models,” _arXiv preprint arXiv:1808.02608_ , 2018.
  * [32] A. Nagrani, J. S. Chung, and A. Zisserman, “VoxCeleb: A Large-Scale Speaker Identification Dataset,” _Proc. Interspeech 2017_ , pp. 2616–2620, 2017.
  * [33] J. S. Chung, A. Nagrani, and A. Zisserman, “VoxCeleb2: Deep Speaker Recognition,” _Proc. Interspeech 2018_ , pp. 1086–1090, 2018.
  * [34] D. Snyder, D. Garcia-Romero, G. Sell, D. Povey, and S. Khudanpur, “X-vectors: Robust DNN embeddings for speaker recognition,” _Submitted to ICASSP_ , 2018.
  * [35] M. McLaren, L. Ferrer, D. Castan, and A. Lawson, “The Speakers in the Wild (SITW) Speaker Recognition Database,” in _Interspeech_ , 2016, pp. 818–822.
  * [36] Y. Qian, X. Chang, and D. Yu, “Single-channel multi-talker speech recognition with permutation invariant training,” _Speech Communication_ , vol. 104, pp. 1–11, 2018.
  * [37] X. Chang, Y. Qian, K. Yu, and S. Watanabe, “End-to-End Monaural Multi-speaker ASR System without Pretraining,” _arXiv preprint arXiv:1811.02062_ , 2018.


# The impact of ultraviolet heating and cooling on the dynamics and
observability of lava planet atmospheres

T. Giang Nguyen,1 Nicolas B. Cowan,2 Raymond T. Pierrehumbert,3 Roxana E.
Lupu,4 and John E. Moores1
1Centre for Research in Earth and Space Sciences, York University, 4700 Keele
St, Toronto, ON M3J 1P3, Canada
2Department of Earth and Planetary Sciences, and Department of Physics, McGill
University, 3550 Rue University, Montréal, QC H3A 2A7, Canada
3Department of Physics, University of Oxford, Oxford, OX1 3PU, United Kingdom
4BAER Institute/NASA Ames Research Center, Moffet Field, CA 94035, USA E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Lava planets have non-global, condensible atmospheres similar to icy bodies
within the solar system. Because they depend on interior dynamics, studying
the atmospheres of lava planets can lead to understanding unique geological
processes driven by their extreme environment. Models of lava planet
atmospheres have thus far focused on either radiative transfer or
hydrodynamics. In this study, we couple the two processes by introducing
ultraviolet and infrared radiation to a turbulent boundary layer model. We
also test the effect of different vertical temperature profiles on atmospheric
dynamics. Results from the model show that UV radiation affects the atmosphere
much more than IR. UV heating and cooling work together to produce a
horizontally isothermal atmosphere away from the sub-stellar point regardless
of the vertical temperature profile. We also find that stronger temperature
inversions induce stronger winds and hence cool the atmosphere. Our simulated
transmission spectra of the bound atmosphere show a strong SiO feature in the
UV that would be challenging to observe in the planet’s transit spectrum due
to the precision required. Our simulated emission spectra are more promising,
with significant SiO spectral features at 4.5 and 9 $\mu$m that can be
observed with the James Webb Space Telescope. Different vertical temperature
profiles produce discernible dayside emission spectra, but not in the way one
would expect.

###### keywords:

instrumentations: detectors, methods: numerical, planets and satellites:
atmospheres

††pubyear: 2021††pagerange: The impact of ultraviolet heating and cooling on
the dynamics and observability of lava planet atmospheres–A

##  1 Introduction

###  1.1 Background

Lava planets are rocky exoplanets that orbit close to their star. Being
tidally-locked into synchronous rotation, these planets have large surface
temperature contrasts between the permanent dayside and nightside. On the
dayside, the stellar flux is intense enough to vaporize rocks, generating a
thin atmosphere of mineral vapour (Schaefer & Fegley, 2009). On the nightside,
the temperature drops and the atmosphere may collapse back onto the surface as
has been predicted for volatiles on M-Earths (Wordsworth, 2015). In order to
qualify as a lava planet, a world must have a density consistent with rocky
bulk composition, and the surface temperature must exceed the solidus for
silicates; examples include CoRoT-7b (Léger et al., 2011), Kepler-10b (Batalha
et al., 2011), and K2-141b (Malavolta et al., 2018).

Dynamics of the partial atmospheres on lava planets are similar to that of Io
(Tsang et al., 2016) or Pluto (Gladstone & Young, 2019) where pressure can
vary significantly in response to huge changes in temperature. As such, icy
bodies in the solar system are counter-intuitively analogous to lava planets
and studying one can help to understand the other. The atmosphere on a lava
planet is in vapour equilibrium with the magma ocean, so atmospheric and
interior dynamics are intertwined (Kite et al., 2016). Therefore, atmospheric
detection and characterization of a lava planet may reveal the bulk
composition of its rocky mantle.

Ito et al. (2015) used radiative transfer simulation to show that a silicon
monoxide atmosphere near the substellar point can be much hotter than the
surface and exhibit a strong inversion (negative lapse rate). Infrared
spectroscopy should reveal dayside SiO emission features. However, such 1D
radiative transfer simulations neglect atmospheric circulation.

In contrast, Nguyen et al. (2020) focused on the hydrodynamical flow of the
silicate atmosphere, with an emphasis on the horizontal advection of heat and
the spatial extent of the atmosphere. We had modelled an optically thin
silicate atmosphere that did not absorb or emit any radiation. This made it
difficult to simulate observables such as emission and transmission spectra.
By construction, the optically thin atmosphere does not absorb stellar
radiation and thus was everywhere colder than the surface. A temperature
inversion as envisioned by Ito et al. (2015) might affect the atmospheric
dynamics. Therefore, the next step is to combine radiative transfer and
atmospheric circulation to better understand lava planets atmospheric dynamics
and more faithfully simulate observations.

###  1.2 Objectives

Castan & Menou (2011) and Nguyen et al. (2020) use the same formulation — they
are insightful for atmospheric dynamics but both assume a transparent
atmosphere. Since SiO vapour absorbs infrared and ultraviolet radiation very
well, we aim to eliminate this approximation by adding radiative transfer to
the evaporation-driven hydrodynamics. Hydrodynamical simulations have also
assumed a dry adiabatic vertical temperature profile, contrary to the 1D
radiative transfer simulations of Ito et al. (2015). We therefore also test a
wider range of vertical temperature profiles to harmonize with the radiative
transfer simulations.

The resulting model couples radiative transfer to atmospheric circulation,
something yet to be done for lava planets. As an ancillary benefit, they also
allow us to more faithfully predict observables. By calculating the
transmission and emission spectra, we simulate JWST observations. We continue
using K2-141b as our poster child as it currently has the highest signal-to-
noise ratio among known lava planets making it the best target for atmospheric
detection (Zieba et al., 2022) and is one of the first planets that JWST will
look at (Dang et al., 2021).

The next section, Section 2, describes the model’s formulation: introduction
of the general turbulent boundary layer (2.1), simplification of the general
equations to 1D (2.2), and the implementation of radiative transfer (2.3).
Section 3 presents results from the model. Section 4 discusses the impact of
UV radiation (4.1.1) and different vertical temperature profiles (4.1.2) on
atmospheric dynamics, as well as presenting simulated transit spectra (4.2.1)
and eclipse spectra (4.2.2). We conclude in Section 5.

##  2 Theory

###  2.1 The turbulent boundary layer

Following the formulation of Ingersoll et al. (1985), the atmosphere consists
of a turbulent boundary layer that behaves similarly to the shallow-water
equations. The surface pressure is assumed to be in fast equilibration with
the magma ocean. As the flow is assumed to be in steady-state, we have the
following general equations describing the conservation of mass, momentum, and
energy respectively:

$\nabla\cdot(\rho hV)=m\ E,$ (1) $\nabla\cdot(\rho
hV^{2})=-\nabla\int_{z}Pdz+\tau,$ (2) $\nabla\cdot\bigg{(}\rho
hV\left(\frac{V^{2}}{2}+C_{p}T\right)\bigg{)}=Q,$ (3)

where $\rho$ is the air density, $h$ is the atmospheric column thickness, $V$
is wind velocity, $m$ is the mass per molecule, $P$ is the pressure, $\tau$ is
the surface drag, $C_{p}$ is the heat capacity, $T$ is the atmospheric
temperature, and $Q$ is the net energy flux for the boundary layer. The
evaporation/outgassing rate $E$, is calculated via:

$E=\frac{P-P_{v}(T_{s})}{m\sqrt{2\pi RT_{s}}},$ (4)

where $R$ is the gas constant, $T_{s}$ is the surface temperature, and $P_{v}$
is the temperature-dependent saturation vapour pressure and can be found in
Miguel et al. (2011).

The surface drag, $\tau$, is calculated by:

$\tau=-\rho Vw,$ (5)

where $w$ is the transfer coefficient defined by Ingersoll et al. (1985) and
is further parameterized by the mean flow velocity, $V_{e}=mE/\rho$, and the
eddy velocity, $V_{d}=V_{*}^{2}/V$. The frictional velocity, $V_{*}$, is
implicitly solved from the experimental formulation of turbulent flow over
smooth flat plates:

$V=2.5V_{*}\log\bigg{(}\frac{9.0\ V_{*}\ H\rho}{2\eta}\bigg{)},$ (6)

where $H$ is the atmospheric scale height and $\eta$ is the dynamic viscosity.
We can now compute the transfer coefficient $w$ for the two cases where there
is net evaporation ($V_{e}>0$) and net deposition ($V_{e}\leq 0$):

$\displaystyle w$ $\displaystyle=\frac{2V_{d}^{2}}{V_{e}+2V_{d}}$
$\displaystyle(V_{e}>0)$ (7a) $\displaystyle w$
$\displaystyle=\frac{V_{e}^{2}-2V_{d}V_{e}+2V_{d}^{2}}{-V_{e}+2V_{d}}$
$\displaystyle(V_{e}\leq 0).$ (7b)

As the flow is assumed to be turbulent, $V$ can be vertically uniform. To
solve the momentum equation, $P$ is vertically integrated with height by using
the pressure and temperature at the top of the boundary layer, $P_{*}$ and
$T_{*}$ respectively:

$\int_{0}^{\infty}Pdz=\int_{0}^{P_{*}}\frac{\Phi}{g}dP=\frac{C_{p}}{g}\int_{0}^{P_{*}}(T_{*}-T)dP,$
(8)

where $\Phi$ is the geopotential. The change in $\Phi$ is the work done as the
vaporized particles ascend from the surface to the top of the boundary layer;
this can be expressed as $C_{p}(T_{*}-T)$. In order to close the problem, we
must assume a vertical temperature profile. For an adiabatic profile, we have
$T=T_{*}(P/P_{*})^{R/C_{p}}$. Substituting this into the integral in Eq. 8
yields:

$\frac{C_{p}}{g}\bigg{(}T_{*}P_{*}-T_{*}P_{*}\bigg{(}\frac{1}{R/C_{p}+1}\bigg{)}\bigg{)}=\frac{C_{p}T_{*}P_{*}\beta}{g},$
(9)

where $\beta=R/(R+C_{p})$ for an adiabatic profile. For an isothermal
atmosphere, on the other hand, $\beta=1$ . Likewise an atmosphere with a
temperature inversion would have $\beta$ greater than 1 (Figure 1 shows
temperature – pressure profiles for different $\beta$). The isothermal and
inverted profiles are more in line with the temperature profile predicted by
Ito et al. (2015). We will test these three different profiles and see their
effects on the dynamics and observability of the K2-141b’s atmosphere.

Figure 1: T-P profiles at the substellar point adopted by Nguyen et al. (2020)
for different $\beta$. The black line represents the adiabatic profile. The
red line represents isothermal and the two blue lines represent profiles with
a temperature inversion.

Because of our defined shallow turbulent layer, the vertical integration of
pressure, temperature, and wind speed allow us to evaluate these state
variables at the layer’s boundary. Vertical integration requires an assumption
of the temperature profile which we parameterized with $\beta$. The
calculations become akin to the shallow-water equations where height is fixed
within the column.

###  2.2 Hydrodynamics

We must now reduce the general equations 1-3 to a more manageable set of
calculations. As K2-141b is tidally locked, we can impose axial symmetry and
the problem reduces to 2 dimensions: height and angular distance from the
substellar point $\theta$. Assuming a shallow turbulent boundary layer further
reduces height dependence. The problem becomes 1D where we solve for pressure,
temperature, and wind velocity with respect to $\theta$ at the top of the
turbulent layer. The general equations describing the conservation of mass,
momentum, and energy become:

$\frac{1}{r\ \sin(\theta)}\frac{d}{d\theta}\bigg{(}\frac{V_{*}\ P_{*}\
\sin(\theta)}{g}\bigg{)}=m\ E,$ (10) $\frac{1}{r\
\sin(\theta)}\frac{d}{d\theta}\bigg{(}\frac{(V_{*}^{2}+\beta\ C_{p}\
T_{*})P_{*}\ \sin(\theta)}{g}\bigg{)}=\frac{\beta\ C_{p}\ T_{*}\ P_{*}}{g\ r\
\tan(\theta)}+\tau,$ (11) $\frac{1}{r\
\sin(\theta)}\frac{d}{d\theta}\bigg{(}\frac{(V_{*}^{2}/2+\ C_{p}\ T_{*})V_{*}\
P_{*}\ \sin(\theta)}{g}\bigg{)}=Q,$ (12)

where K2-141b’s radius, $r$, and surface gravity, $g$, are $9.62\times 10^{6}$
m and $21.8$ m s-2 respectively.

A limitation of the model is that the atmosphere’s composition must be pure.
Due to its volatility, sodium has been the prime outgassing candidate for
previous lava planet modelling papers (Castan & Menou, 2011; Kite et al.,
2016). However, Na is relatively scarce in the crust, and risks being lost to
space (Ito & Ikoma, 2021) or cold trapped on the nightside (Nguyen et al.,
2020), so we argue instead for a silicate atmosphere due to the abundance of
SiO2 in the crust. Both SiO2 and SiO can be evaporated into the atmosphere but
the saturated vapour pressure of SiO2 is several orders of magnitude lower
(Schaefer & Fegley, 2009). Therefore, the atmosphere is likely dominated by
SiO and its relevant parameters are: $m=7.32\times 10^{-26}$ kg molecule-1,
$C_{p}=851$ J K-1 kg-1, $R=188$ J K-1 kg-1.

The system of ordinary differential equations, Eq. 10-12, can be solved
numerically as a boundary value problem. To find a solution, we use the
shooting method where we test different sublimation rates at the substellar
point to see which value leads to steady-state flow. When including radiative
transfer, as described below, we must also test initial energy flux because
evaporation is no longer the sole driving force of the flow.

###  2.3 Radiative transfer

To implement radiative transfer in our model, we need the wavelength-dependent
absorptivity of SiO and the stellar spectrum, both shown in Fig. 2. SiO
absorbs strongly in the ultraviolet and relatively weakly in the infrared.
Since SiO has negligible optical absorption, we only consider the IR and UV
bands for heating and cooling of the atmosphere.

Figure 2: _Top_ : absorption cross-section of SiO at 1 kPa and 2800 K. The UV
data are from Kurucz (1992) while the remainder are from Barton et al. (2013).
Broadening parameters are approximated based on HITRAN data for similar
diatomic molecules, assuming an H2/He background atmosphere (Gordon et al.,
2017); although this assumption is hard to justify, this is currently the best
source of SiO absorption cross-section at high temperatures. _Bottom_ :
stellar spectrum of HD 5512 (France et al., 2016), a close analogue of K2-141.

Implementing radiative transfer starts with determining the absorptivity, or
equivalently the emissivity, $\epsilon$, in each waveband:

$\epsilon=1-e^{-\tau_{d}},$ (13)

where $\tau_{d}$ is the optical depth. This depth is determined using the
surface pressure:

$\tau_{d}=xP_{*}/mg,$ (14)

where $x$ is the absorption cross-section of SiO in the IR or UV: $x_{\rm
IR}=10^{-29}$ m2/mole and $x_{\rm UV}=10^{-22}$ m2/mole.

With emissivity and absorptivity accounted for, we can now define the
atmospheric heating $Q$ as:

$Q=Q_{\rm sens}+F_{\rm IR}+F_{\rm UV}+F_{\rm surf}-2F_{\rm RC},$ (15)

where $Q_{\rm sens}$ is the sensible heating as described by Ingersoll et al.
(1985), $F_{\rm IR}$ is the atmospheric absorption of stellar IR, $F_{\rm UV}$
is the atmospheric absorption of stellar UV, $F_{\rm surf}$ is the blackbody
radiation emitted by the surface and subsequently absorbed by the atmosphere,
and $F_{\rm RC}$ is the radiative cooling of the atmosphere. The radiative
heating terms for the atmosphere are:

$F_{\rm IR}=\epsilon_{\rm IR}\ C_{\rm IR}\ F_{*},$ (16) $F_{\rm
UV}=\epsilon_{\rm UV}\ C_{\rm UV}\ F_{*},$ (17) $F_{\rm surf}=\epsilon_{\rm
IR}\sigma T_{s}^{4}.$ (18)

In the equations above, $F_{*}$ is the incident stellar flux, which is a
function of $\theta$ calculated following Kopal (1954) to account for the non-
negligible angular size of the star as seen from the planet. For the stellar
spectrum, we use the star HD 5512 as it is a close analogue of K2-141; $C_{\rm
IR}=0.532$ and $C_{\rm UV}=0.009$ are the fraction of stellar flux in the IR
and UV range, respectively (France et al., 2016). We ignore stellar UV
absorption via photo-dissociation of SiO because, as we show below,
absorptivity and emissivity in the UV is $\sim$1 for most the dayside. Since
nearly all of the stellar UV has been absorbed, any processes that can further
absorb UV are negligible. Furthermore, significant photo-dissociation of SiO
occurs at wavelengths less than 0.14 $\mu$m (Jolicard et al., 1997) and
stellar flux in this range is negligible.

For radiative cooling, only considering IR emission would yield $F_{\rm
RC}=\epsilon_{\rm IR}\sigma T_{*}^{4}$. However, since the absorption cross-
section of the atmosphere is several orders of magnitude larger in the UV than
the IR, UV emission becomes non-negligible. Therefore we split the
$\epsilon\sigma T_{*}^{4}$ term using definite integrals of the Planck
function, $B$, over IR and UV wavelengths, $\lambda$, and the equation for
$F_{\rm RC}$ becomes:

$F_{\rm RC}(T)=\pi\bigg{(}\epsilon_{\rm IR}\int_{\lambda_{\rm
IR}}B(\lambda,T)d\lambda+\epsilon_{\rm UV}\int_{\lambda_{\rm
UV}}B(\lambda,T)d\lambda\bigg{)}.$ (19)

A computationally efficient way to evaluate the Planck integral can be found
in the Appendix. The surface albedo is assumed to be 0 and this is justified
as quenched glass and liquid lava surfaces are expected to have very low
albedo (Essack et al., 2020). We account for the greenhouse effect by
radiatively coupling the atmospheric and surface temperatures:

$\sigma T_{s}^{4}=F_{*}-F_{\rm IR}-F_{\rm UV}+F_{\rm RC}(T).$ (20)

##  3 Results

Starting with an adiabatic profile ($\beta=0.18$), we first only include IR
absorption and emission. The results did not differ much from a transparent
atmosphere aside from having slightly stronger winds and warmer atmospheric
temperatures. Atmospheric temperatures are still always cooler than the
surface (cf. dotted and dashed lines in Fig. 3). Pressure almost never changes
regardless of the radiative schemes or temperature profiles used as it cannot
deviate far from the local saturation vapour pressure.

Figure 3: Atmospheric pressure, wind velocity and temperature from 1D
hydrodynamic + radiative transfer simulations of K2-141b assuming an adiabatic
vertical temperature-pressure profile. In each panel the dotted line
represents the transparent atmosphere (no radiative transfer, as in Nguyen et
al., 2020), the dashed line represents the simulation with only IR radiative
transfer, and the solid line is the simulation incorporating both IR and UV
radiative transfer. The IR-only radiative scheme is similar to the transparent
case but the scheme with UV radiation produces stronger winds and a much
hotter atmosphere. It can be seen that the temperature for the IR-only case
rises significantly at high $\theta$. This is caused by viscosity which slow
winds and raises the temperature, but is only effective at low temperatures.
However, the transparent atmosphere is just as cold but does not heat up. This
is because it has significantly less pressure than its IR-only counterpart and
there is not enough atmospheric mass for friction to take effect.

When we add UV radiation in addition to IR, the resulting winds are stronger
and the atmospheric temperature exceeds that of the surface for
$\theta>32^{\circ}$ (solid lines in the bottom panel of Fig. 3). The reason
for the high atmospheric temperature is the high optical depth in the UV
compared to that in the IR. As shown in the middle panel of Fig. 4, UV
absorptivity hovers around unity while the IR opacity, being very sensitive to
pressure, drops exponentially with $\theta$.

The bottom panel of Fig. 4 shows that IR absorption from stellar flux and
surface blackbody radiation start out as the dominant heating terms, balanced
by IR radiative cooling. As pressure drops, however, IR radiation becomes
negligible and the UV radiative terms come to dominate the energy balance.

Figure 4: Radiative properties of the adiabatic atmosphere. _Top_ : vertical
optical depth of the SiO atmosphere. _Middle_ : vertically integrated
absorptivity (equivalently emissivity). _Bottom_ : radiative fluxes absorbed
and emitted by the atmosphere. The emissivity in the UV is $\sim$1 for most of
the dayside while emissivity in the IR drops exponentially away from the
substellar point. While UV heating and cooling are relative small near the
sub-stellar point, they dominate the radiation budget at $\theta>50^{\circ}$.

Finally, we incorporate IR and UV radiation to two other vertical temperature
profiles: isothermal ($\beta=1$) and inverted ($\beta=1.2$). Pressure, wind
speed, and temperature for the three vertical profiles are shown in Fig. 5.
Higher $\beta$ produce stronger winds and a cooler atmosphere; in Section
4.1.2, we offer a physical explanation of this phenomenon. Near the day-night
terminator, however, atmospheric temperatures are the same regardless of
$\beta$; we discuss this in Section 4.1.1.

Figure 5: Results of hydrodynamic simulations including IR and UV radiative
transfer for different vertical temperature profiles. _Top_ : atmospheric
pressure. _Middle_ : wind velocity. _Bottom_ : atmospheric temperature. The
black line is the adiabatic profile, red is isothermal, and blue is the
inverted temperature profile. Larger $\beta$ increase wind speed and decrease
temperature. However, the temperature for all three profiles converge at
$\theta>80^{\circ}$.

##  4 Discussion

###  4.1 Impact on dynamics

####  4.1.1 UV radiation

Our simulations show that the exchange of IR radiation between surface and
atmosphere does not change the dynamics significantly despite the large
fluxes. This is because IR heating and cooling are well balanced, leading to a
near zero net radiation budget. However, IR radiation becomes negligible at
low pressure making UV stellar flux the dominant heating term, balanced by UV
radiative cooling.

The importance of UV cooling in a planetary atmosphere is surprising but
understandable. Fig. 6 shows the IR and UV cooling at a fixed pressure. For
uniform emissivity at all wavelengths, the blackbody radiation in the IR
accounts for almost all of the radiative cooling. However, because emissivity
in the IR is much lower than in the UV, UV emission can overtake the IR with
only a slight increase in temperature.

Figure 6: _Top_ : atmospheric radiative cooling terms featured in Eq. 19. The
black line is the total blackbody radiative flux across all wavelengths, red
is the IR contribution, and blue is the UV contribution. At constant
emissivity, IR flux would be order of magnitude greater than UV flux. _Bottom_
: emitted blackbody radiation at a fixed pressure. The dashed black line is
the total emitted flux, the solid black line is what the emitted flux would be
if the total emission is dependent only on IR emissivity, the red line is the
individual emitted flux in the IR and blue in the UV. Note that UV emission
becomes non-negligible when atmospheric temperatures exceed 2500 K.

The temperature sensitivity of UV emission explains the horizontal uniformity
in temperature, especially near the terminator ($\theta=90^{\circ}$). As
pressure drops, UV stellar flux is the only heating term and IR cooling is
ineffective. The atmosphere will get hotter until UV radiative cooling takes
over. Therefore, temperature will hover around 2600 K where UV heating and
cooling are in equilibrium, regardless of the assumed vertical profile (see
bottom panel of Fig. 5).

####  4.1.2 Vertical temperature profile

As described in section 2.1, the model atmosphere’s vertical temperature
profile can be made isothermal by setting $\beta=1$ or inverted for $\beta>1$.
This stems from the vertical integration of the pressure which only affects
the momentum equation (Eq. 11). The higher $\beta$ value implies a much
stronger pressure force, the term $\beta C_{p}T_{*}P_{*}/[gr\tan(\theta)]$ in
Eq. 11. This is why winds are stronger for isothermal and inverted temperature
profiles as seen in Fig. 5. We plot the temperature profiles at two different
locations in Fig 7.

Figure 7: Temperature profiles at $\theta=30^{\circ}$ (left) and
$\theta=45^{\circ}$ (right). While the adiabatic profile stays largely the
same at both location (cf. Bottom panel of Fig. 5), the isothermal and
inverted temperature profiles are cooler at $\theta=45^{\circ}$.

Stronger winds also induce more evaporation and condensation and this can be
inferred from the mass equation (Eq. 10): a larger $V_{*}$ must lead to a
larger $E$. Therefore, a strong vertical temperature inversion has a greater
evaporation and condensation rate as seen in Fig. 8.

Figure 8: Evaporation rate for different vertical temperature profiles.
Because decreasing the lapse rate results in stronger winds, evaporation and
condensation are greater for the isothermal and inverted temperature profiles
than the adiabatic profile.

Wind speed and atmospheric temperature are coupled within the energy equation
(Eq. 12), hence the increase in wind speed greatly affects temperature. The
term $(V^{2}/2+C_{p}T_{*})$ is balanced by $Q$ which is mostly dependent on
the incoming stellar flux. As stellar flux does not change regardless of the
vertical temperature profile, a lower temperature is required to balance out
the stronger wind speeds.

Radiative transfer models that do not account for horizontal flow predict an
inverted temperature profile where the overall atmospheric temperature is
hotter than the surface (Ito et al., 2015). The hydrodynamics we present here
suggests otherwise: the stronger the inversion, the faster the horizontal flow
and the cooler the atmosphere. However, because the lapse rate is set up to be
the same everywhere ($\beta$ is constant with respect to $\theta$), the
hydrodynamical effects may be overestimated.

###  4.2 Impact on observations

####  4.2.1 Transmission spectroscopy

We can now use our model atmospheres to simulate observations, the first of
which is a transit spectrum. This is done by calculating the transit depth,
$D$, the fraction of stellar flux blocked as the planet passes in front of its
star. This is evaluated by:

$D(\lambda)=\bigg{(}\frac{r}{r_{*}}\bigg{)}^{2}+\frac{2}{r_{*}^{2}}\int_{r}^{r_{*}}b(1-e^{-\tau_{d}(\lambda,b)})db,$
(21)

where $r$ and $r_{*}$ are the planetary and stellar radii respectively, $b$ is
the projected distance to the planet’s centre, and $\tau_{d}$ is the optical
depth along the cord. Optical depth is calculated by:

$\tau_{d}(b)=\frac{xP_{*}}{mg}\sqrt{\frac{2\pi r}{H}},$ (22)

where $H$ is the atmospheric scale height. In the middle of the transit, the
relevant surface pressure $P_{*}$ and $H$ must be evaluated at
$\theta=90^{\circ}$, but the extreme geometry of lava planets means that the a
wide range of $\theta$, surface pressures, and scale heights are probed
throughout transit, as shown in Fig. 9. Because the boundary layer is
geometrically thin compared to scale height, $P_{*}$ can be approximated as
the surface pressure.

Due to the UV radiative equilibrium, atmospheric temperatures converge to the
same temperature at the terminator regardless of the vertical temperature
profile. These temperatures are much hotter than that of the transparent
atmosphere described by Nguyen et al. (2020), increasing the atmospheric scale
height and providing a stronger signal for transmission spectroscopy of lava
planets.

Evaluating $\tau_{d}$ and completing the integration in Eq. 21 produces the
transit spectra shown in Fig. 10. The biggest peak is at $\sim$0.25 $~{}\mu$m
and the second biggest at $\sim$9$~{}\mu$m, both corresponding to SiO
absorption features since our model atmosphere is entirely composed of that
gas. By virtue of having slightly higher atmospheric pressure and therefore
greater optical depth at the terminator (see Fig. 5), the isothermal and
inverted atmospheres have larger transit spectrum features than the adiabatic
case.

The spectral feature in the UV is the strongest but likely requires precision
beyond the capabilities of the Hubble Space Telescope. But we only simulate
the bound atmosphere of the planet—it is possible that UV observations with
Hubble would be sensitive to its exosphere. Transit spectral features in the
range of JWST NIRSpec and Spitzer IRAC are negligible while features in the
range of JWST MIRI are significant. We use the PandExo tool from Batalha et
al. (2017) to simulate JWST measurement uncertainties for 16 transits and
conclude that neither NIRSpec or MIRI have the precision to identify any
spectral features induced by the atmosphere in transit.

Figure 9: Scale height as observed at ingress (left) and transit (right).
Assuming an impact parameter of 0, we plot the expected scale height around
the planet’s limb. At ingress, the viewing geometry exposes much of the
airless nightside (left side of the limb) and we can only see the dayside
atmosphere on the right-hand side. At transit, limb measurements only probe
where $\theta=90^{\circ}$ and scale height is constant due to the imposed
axial symmetry. Incorporating radiative transfer in our dynamical model
produces a much hotter and hence more extended atmosphere near the terminator,
leading to a near-constant transmission spectrum throughout the transit.
Figure 10: Expected transit spectrum of K2-141b for three possible vertical
temperature profiles. Dashed lines show the bandpasses for HST/COS (Green et
al., 2011), JWST/NIRSpec (Birkmann et al., 2016), Spitzer/IRAC channel 2, and
JWST/MIRI (Wells et al., 2015). The UV spectral feature is in-band for Hubble
but is likely too small to be detectable. Likewise, the IR features around 4.5
and 9 $\mu$m, while in-band for Spitzer or JWST, are too faint to be
disentangled from noise. Although cases with different T-P profiles have
roughly the same scale height at $\theta=90^{\circ}$, the surface pressure is
slightly larger for profiles with a negative lapse rate. This explains the
slight differences in transit depth between the T-P profiles.

####  4.2.2 Emission spectroscopy

To calculate the emission spectrum of K2-141b at eclipse, we first calculate
the top-of-atmosphere outgoing irradiance, which is a combination of
atmospheric and planetary emission. While the greenhouse effect barely changes
the surface temperature, it is significant enough to affect the planet’s
emission spectrum. Emission from the atmosphere is calculated as
$\epsilon_{\lambda}B(\lambda,T_{*})$ while the surface contribution is
calculated as $(1-\epsilon_{\lambda})B(\lambda,T_{s})$. The top two panels of
Fig. 11 show the atmospheric and surface irradiance at 4.5 $\mu$m, where there
is a strong SiO spectral feature. The bottom panel shows the total top-of-
atmosphere outgoing irradiance.

Figure 11: Irradiance of K2-141b at 4.5 $\mu$m i.e, at a wavelength where the
atmosphere is optically thick over much of the dayside. _Top_ : modeled
outgoing atmospheric irradiance. _Middle_ : outgoing planetary irradiance.
_Bottom_ : total outgoing irradiance at the top of the atmosphere. The
adiabatic case has stronger emission at this wavelength while the other two
cases have a weaker emission than the transparent atmosphere.

The emission spectrum is a function of orbital phase because different regions
of the planet rotate in and out of view. Therefore, outgoing irradiance with
respect to $\theta$ is translated to a spherical coordinates with co-latitude
($\alpha$) and longitude ($\phi$). We then integrate the irradiance over the
dayside hemisphere following Cowan & Agol (2008) to find the emission,
$I_{e}$:

$I_{e}(\lambda,\zeta)=\int_{0}^{\pi}\int_{\zeta-\pi/2}^{\zeta+\pi/2}I(\alpha,\phi,\lambda)\sin^{2}(\alpha)\cos(\phi-\zeta)d\phi
d\alpha,$ (23)

where $\zeta$ is the position of the planet in its orbit: $\zeta=0$ at
superior conjunction and $\zeta=\pi$ at inferior conjunction. Keeping
$\lambda$ fixed at 4.5 $\mu$m and calculating $I_{e}(\zeta)$ yields the phase
curve shown in Fig. 12. Because the dayside emission from the isothermal and
inverted temperature profiles are smaller than the adiabatic case, their phase
amplitudes are likewise smaller. Eclipse occurs when $\zeta=0$, and
calculating $I_{e}$ with respect to $\lambda$ at this point leads to the
eclipse spectrum plotted in Fig. 13.

Paradoxically, the adiabatic profile results in SiO emission features while
the isothermal and inverted profiles result in absorption features. SiO
spectral features probe the atmosphere rather than the surface over much of
the dayside. The adiabatic profile produces an atmosphere hotter than the
surface and therefore its spectrum shows emission features. Meanwhile, the
isothermal and inverted temperature profiles lead to cooler atmospheres
(bottom panel of Fig. 5) and hence their spectra show absorption features.

We simulate JWST measurement uncertainties for 16 eclipses using PandExo
(Batalha et al., 2017), plotted in the bottom panel of Fig. 13. Spectral
features at 4.5 $\mu$m can be detected with NIRSpec. Spectral features at
longer wavelengths are difficult to see with MIRI but observing more eclipses
or phase curves may reduce the errors enough to make out the features.

Figure 12: Simulated 4.5 $\mu$m phase curves from the three different vertical
temperature structures. Note that eclipses have been omitted at phases 0 and
1; transit is omitted at phase 0.5. The adiabatic temperature profile produced
an atmosphere with stronger dayside emission and therefore a larger phase
amplitude. Figure 13: _Top_ : simulated eclipse spectrum of K2-141b for three
different assumed vertical temperature pressure profiles. The black dot and
its bars show K2-141b’s occultation depth as measured by Kepler (Malavolta et
al., 2018), as well as Kepler’s bandwidth (Van Cleve & Caldwell, 2016). The
cyan dot and its bars likewise show the depth as measured by Spitzer (Zieba et
al., 2022). The adiabatic case exhibits SiO emission features while the
vertically isothermal and inverted models have absorption features. This
counter-intuitive behaviour originates from the relation between vertical
temperature profile, wind speed, and overall atmospheric temperature with
respect to the surface. _Bottom_ : simulated JWST observations at eclipse. We
take the spectrum from the isothermal profile case and simulate NIRSpec and
MIRI measurements for 16 eclipses using the tool PandExo (Batalha et al.,
2017). Results for NIRSpec are promising and can detect the SiO absorption
feature at 4–5 $\mu$m.

##  5 Conclusions

We have improved upon the 1D hydrodynamic model of Nguyen et al. (2020) by
implementing radiative transfer. We applied this new model to the lava planet
K2-141b and tested the effects of different vertical temperature profiles. Our
analysis leads to two key insights.

Our first insight is that the atmosphere far from the sub-stellar point will
be hotter than the surface. Since UV heating dominates over IR when pressure
is low, atmospheric temperatures will rise to $\sim$2600K, at which point UV
radiative cooling takes effect. As a result, atmospheric temperatures should
be horizontally uniform near the terminator ($\theta\approx 90^{\circ}$),
regardless of vertical temperature profile.

The second insight is that a vertical temperature inversion produces a cooler
atmosphere. A strong temperature inversion also leads to a greater horizontal
pressure gradient force and hence greater wind speeds. To conserve energy,
thermal energy in the atmosphere is reduced as a reaction to the greater
kinetic energy. This result is particularly important as it shows how
hydrodynamics can lead to a significantly different atmosphere when compared
to models that focus on 1D radiative balance (e.g., Ito et al., 2015).

We also simulated the transit spectrum of K2-141b (Fig. 10) and found that UV
radiation leads to a much bigger scale height at the terminator than what was
previously predicted (Nguyen et al., 2020). However, the surface pressure is
very small and so features in the IR are too subtle for Spitzer or JWST to
detect, while the larger UV SiO feature likely have to wait for a more
sensitive spectrograph than is currently available on Hubble.

Our simulated eclipse spectra of K2-141b (Fig. 13) are more promising,
however. Because the different vertical temperature profiles lead to different
dayside atmospheric temperatures, emission spectra ought to be distinguishable
with Spitzer and JWST. The adiabatic profile has a warmer atmosphere so its
spectrum has significant SiO emission features. The isothermal and inverted
vertical temperature profiles, on the other hand, produce atmospheres cooler
than the surface and hence SiO absorption features. Emission spectra are not
only useful for detecting the atmosphere, but can also help to infer more
nuanced properties such as vertical temperature profile.

## Acknowledgements

This work was made possible by the Natural Science and Engineering Research
Council (NSERC) of Canada’s Collaborative Research and Training Experience
Program (CREATE) for Technology for Exo-Planetary Sciences (TEPS). This work
was also made possible by the NASA XRP grant NNX17AC02G. This work has made
use of the MUSCLES Treasury Survey High-Level Science Products. NBC
acknowledges support from the McGill Space Institute (MSI) and l’Institut de
recherche sur les exoplanètes (iREx).

## Data Availability

Data is available upon request.

## References

  * Barton et al. (2013) Barton E. J., Yurchenko S. N., Tennyson J., 2013, Monthly Notices of the Royal Astronomical Society, 434, 1469
  * Batalha et al. (2011) Batalha N. M., et al., 2011, The Astrophysical Journal, 729, 27
  * Batalha et al. (2017) Batalha N. E., et al., 2017, Publications of the Astronomical Society of the Pacific, 129, 064501
  * Birkmann et al. (2016) Birkmann S. M., et al., 2016, in Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave. p. 99040B
  * Castan & Menou (2011) Castan T., Menou K., 2011, The Astrophysical Journal Letters, 743, L36
  * Cowan & Agol (2008) Cowan N. B., Agol E., 2008, ApJ, 678, L129
  * Dang et al. (2021) Dang L., et al., 2021, JWST Proposal. Cycle 1, p. 2347
  * Essack et al. (2020) Essack Z., Seager S., Pajusalu M., 2020, The Astrophysical Journal, 898, 160
  * France et al. (2016) France K., et al., 2016, The Astrophysical Journal, 820, 89
  * Gladstone & Young (2019) Gladstone G. R., Young L. A., 2019, Annual Review of Earth and Planetary Sciences, 47, 119
  * Gordon et al. (2017) Gordon I. E., et al., 2017, Journal of Quantitative Spectroscopy and Radiative Transfer, 203, 3
  * Green et al. (2011) Green J. C., et al., 2011, The Astrophysical Journal, 744, 60
  * Ingersoll et al. (1985) Ingersoll A. P., Summers M. E., Schlipf S. G., 1985, Icarus, 64, 375
  * Ito & Ikoma (2021) Ito Y., Ikoma M., 2021, Monthly Notices of the Royal Astronomical Society, 502, 750
  * Ito et al. (2015) Ito Y., Ikoma M., Kawahara H., Nagahara H., Kawashima Y., Nakamoto T., 2015, The Astrophysical Journal, 801, 144
  * Jolicard et al. (1997) Jolicard G., Zucconi J.-M., Drira I., Spielfieldel A., Feautrier N., 1997, The Journal of chemical physics, 106, 10105
  * Kite et al. (2016) Kite E. S., Fegley Jr B., Schaefer L., Gaidos E., 2016, The Astrophysical Journal, 828, 80
  * Kopal (1954) Kopal Z., 1954, Monthly Notices of the Royal Astronomical Society, 114, 101
  * Kurucz (1992) Kurucz R. L., 1992, Revista Mexicana de Astronomia y Astrofisica, vol. 23, 23
  * Léger et al. (2011) Léger A., et al., 2011, Icarus, 213, 1
  * Malavolta et al. (2018) Malavolta L., et al., 2018, The Astronomical Journal, 155, 107
  * Miguel et al. (2011) Miguel Y., Kaltenegger L., Fegley B., Schaefer L., 2011, The Astrophysical Journal Letters, 742, L19
  * Nguyen et al. (2020) Nguyen T. G., Cowan N. B., Banerjee A., Moores J. E., 2020, Monthly Notices of the Royal Astronomical Society, 499, 4605
  * Schaefer & Fegley (2009) Schaefer L., Fegley B., 2009, The Astrophysical Journal Letters, 703, L113
  * Tsang et al. (2016) Tsang C. C., Spencer J. R., Lellouch E., Lopez-Valverde M. A., Richter M. J., 2016, Journal of Geophysical Research: Planets, 121, 1400
  * Van Cleve & Caldwell (2016) Van Cleve J. E., Caldwell D. A., 2016, Kepler Science Document KSCI-19033-002, p. 1
  * Wells et al. (2015) Wells M., et al., 2015, Publications of the Astronomical Society of the Pacific, 127, 646
  * Widger Jr & Woodall (1976) Widger Jr W., Woodall M., 1976, Bulletin of the American Meteorological Society, 57, 1217
  * Wordsworth (2015) Wordsworth R., 2015, The Astrophysical Journal, 806, 180
  * Zieba et al. (2022) Zieba S., et al., 2022, arXiv preprint arXiv:2203.00370

##  Appendix A Definite integral of Planck function

Evaluating the radiative cooling requires a definite integral of the Planck
function (Eq. 19). Doing the integral through finite differences can be
computationally expensive but Widger Jr & Woodall (1976) provide a better
alternative. We start with the definite integral of the Planck function from a
wavelength $\lambda_{1}$ to infinity:

$\int_{\lambda_{1}}^{\infty}\frac{2hc}{\lambda^{5}}\frac{1}{e^{(}hc/\lambda
kT)-1}=-\frac{2\pi k^{4}T^{4}}{h^{3}c^{2}}\sum_{n}Ge^{-nJ},$ (24)

where $G=J^{3}/n+3J^{2}/n^{2}+6J/n^{3}+6/n^{4}$ and $J=hc/(\lambda kT)$.
Therefore, the definite integral between two wavelengths is:

$\int_{\lambda_{1}}^{\lambda_{2}}(...)=\int_{\lambda_{1}}^{\infty}(...)-\int_{\lambda_{2}}^{\infty}(...),$
(25)

and each term is calculated via the summation provided.