Zhuangqun Huang*, Gil Keren*, Ziran Jiang*, Shashank Jain*††thanks: *Equal
contribution.,
David Goss-Grubbs, Nelson Cheng, Farnaz Abtahi, Duc Le, David Zhang,
Antony D'Avirro, Ethan Campbell-Taylor, Jessie Salas, Irina-Elena Veliche, Xi
Chen

# Text Generation with Speech Synthesis for ASR Data Augmentation

###### Abstract

Aiming at reducing the reliance on expensive human annotations, data synthesis
for Automatic Speech Recognition (ASR) has remained an active area of
research. While prior work mainly focuses on synthetic speech generation for
ASR data augmentation, its combination with text generation methods is
considerably less explored. In this work, we explore text augmentation for ASR
using large-scale pre-trained neural networks, and systematically compare
those to traditional text augmentation methods. The generated synthetic texts
are then converted to synthetic speech using a text-to-speech (TTS) system and
added to the ASR training data. In experiments conducted on three datasets, we
find that neural models achieve 9%-15% relative WER improvement and outperform
traditional methods. We conclude that text augmentation, particularly through
modern neural approaches, is a viable tool for improving the accuracy of ASR
systems.

Index Terms: Data augmentation, text synthesis, TTS, ASR

##  1 Introduction

Training an end-to-end automatic speech recognition (ASR) model requires
thousands of hours of transcribed speech data. However, collecting human
annotated speech data is challenging and expensive. Data augmentation using
either labelled or unlabelled data [1, 2] has been used to alleviate this data
scarcity problem. One promising approach is speech data synthesis, which
recently contributed to significant progress in domain adaptation [3],
medication names recognition [4], accurate numeric sequences transcription
[5], low-resource languages [6], etc. Most research in this areas has focused
on text-to-speech (TTS) generation [7, 8, 9, 10] and audio augmentation [11,
12, 13], and their effects on the resulting ASR accuracy. However, the quality
of the texts used for generating synthetic speech data aiming at improving ASR
accuracy has not been sufficiently studied. Text augmentation (TA) can
potentially greatly affect the downstream ASR accuracy, as one may use it to
further diversify and enrich the synthetic speech data. Furthermore, TA
methods may also allow quickly setting up an ASR system for a new domain,
since a text corpus can be more economically collected, and may be sufficient
to reach a target accuracy if effectively augmented.

The existing literature about TA for improving ASR accuracy is limited.
Rosenberg et. al. use a non-neural MaxEnt language model (LM) [14] for TA and
observes a small improvement in word error rate (WER) [9]. Zevallos et. al.
shows a good improvement from the non-augmented dataset [6]. In [6], synthetic
utterances were created by first replacing entities like noun etc. with a
delexicalized form which are fed to a seq2seq model to create diverse
delexicalized paraphrases. The delexicalized form is replaced by filling in
with values found from training data. These methods limit the power of TA for
improved ASR and should be replaced with modern, powerful ones pre-trained on
large text corpora. In addition, other than LibriSpeech, those works only
explore domains where small datasets have simple and short utterances.
Furthermore, the reported improvements are convoluted with contributions from
other components, like careful TTS optimization with in-domain data.

In this work, we improve over the existing literature by performing TA using a
large modern neural network, which benefits from a large-scale corpus for pre-
training, for the purpose of ASR data augmentation. In addition, we provide a
first of its kind rigorous comparison of modern neural network and rule-based
models for the purpose of TA for ASR data augmentation. Specifically, we start
from a text corpus, and compare the impact of different TA techniques on
different domains with distinct utterance styles. We use both a rule-based
non-neural approach (NLX), and two other neural methods. Our neural model
based TA methods utilize BART, a large pre-trained model [15] which has proven
useful in large scale synthetic text generation with minimal human effort
[16]. We utilize mask filling to create synthetic utterances where it can find
similar words from the large text corpus it was pre-trained on, unlike [6]
which finds words only from the smaller seed data training corpus. These TA
methods are compared on the open domain LibriSpeech, a Photo Search domain
with short, structured utterances, and a Question and Answering (Q&A) domain
with relatively free-style dictations. To better compare the TA effects, in-
domain data from both the Photo Search and Q&A domains are synthetically
generated, and no domain-specific optimization of the TTS systems is
conducted. Our key findings and contributions are as follows

  * •

TA considerably improves ASR accuracy for open domains or domains with diverse
utterances with 10%-16% relative WER improvements.

  * •

Modern neural based TA is superior to the rule-based approach in terms of
downstream WER gains for all the three domains investigated. As this finding
does not appear in previous literature, it suggests that neural based TA
together with a TTS system is an important viable tool for improving ASR
accuracy.

  * •

For domain-specific training, it is possible to start from a small corpus of
in-domain text, and modern TA technique should be considered as a standard
practice in speech data synthesis for improving ASR systems.

  * •

Even using only the seed text data for TTS data augmentation with no TA
results in significant WER improvements, e.g., relatively $\sim$5% for open
domains and $\sim$15% for domains with simple, structured utterances.

  * •

The model improvement is insensitive to the sizes of augmented data after 8x
augmentation, unless reaching a $>$ 100x augmentation factor.

  * •

Combining synthetic data generated from multiple TA methods further improves
WER gains for the Q&A domain, e.g. 11.88% vs. 15.68% relatively from the
baseline.

##  2 Data

###  2.1 Datasets

We experiment with the LibriSpeech dataset [17] consisting of 960 hours of
transcribed training utterances. In this work, both the LibriSpeech test-clean
and test-other splits in the evaluation dataset are used for model testing. We
also explore two in-house datasets, namely for a Photo Search and a Q&A ASR
domains, where the baseline corpus includes 25k hours of speech from the voice
assistant domains and another 125k hours from videos publicly shared by our
Facebook users. For the Q&A domain, we use two internal in-domain test sets
that we name "in-domain-1" and "in-domain-2", each containing 6k utterances,
originating from live traffic data contributed by internal Meta employees and
from a third-party vendor, respectively.

The in-domain test set for Photo Search consists of 16.2k utterances from
crowdsourcing. For this domain, the instructions in data collection involve
issuing utterances to search and retrieve images or query information about
them, using either an actual device or a simulator. The collected data is then
annotated and marked by humans as success or failure based on the outcome of
the query.

###  2.2 Text Corpus for Augmentation

For text generation, we start from a seed corpus that serves as the input to
each of the generation methods. For LibriSpeech, the seed corpus contains all
the reference text from 281k utterances in the training split. For the Q&A
domain, we use 100k sentences that were gathered by paraphrasing Q&A user
requests from Facebook groups. The seed corpus for the Photo Search domain
contains only 3350 utterances collected using human annotators. The Photo
Search domain has well-structured short utterances, e.g. “show me photos from
new year’s”, while the Q&A domain contains free-style dictations, e.g., “could
I get some suggestions on ingredients that would go well with white rice to
give it a pop of taste”, etc.

##  3 Method

The general pipeline starts from a seed text corpus that is used by a text
augmentation method to emit a new corpus of augmented texts. A TTS system then
takes the augmented texts and generates synthetic speech utterances. We use a
recurrent neural network transducer (RNN-T) model for ASR described in detail
in Section 3.4. During RNN-T training, both the human annotated data and the
augmented TTS data are used, with a hyperparameter controlling their mixing
ratio. We train the RNN-T model from scratch for all but the Photo Search
domain, where the seed data is small and experiments are done by fine tuning
an existing model trained on the data described in Section 2.1.

###  3.1 Text Augmentation

We compare and contrast both rule-based and model-based approaches to text
generation. Early natural language generation (NLG) systems [18] utilized
content selection, macro / micro planning and surface realization. These were
used to create dialog data, whereas our goal is to generate plain text
utterances.

####  3.1.1 NLX

Our rule-based method, called NLX, is a feature-enhanced context-free rule-
based system developed in-house. For NLX, a human developer analyzes the seed
data to produce syntax rules and uses various tools to mine lexical items such
as named entities. The grammar is then used to synthesize a set of utterances.

####  3.1.2 BART

Our model based approach uses BART [15], a neural network-based model. For
BART, we follow the generation scheme from [16] where the vanilla BART model
is first fine-tuned on masked sequences from seed training data which helps
with task adaptation similar to the approach in [19]. During generation we run
a mask filling approach to generate synthetic texts. We use two kinds of
masking strategies: 1) Random Masking (BART R.M.), where we insert masks to
seed data randomly. 2) Custom Masking (BART C.M.), where we use Part of Speech
tags to identify candidate words for masking like nouns and verbs, and replace
them one at a time with masks.

Table 1: Samples of synthetic texts generated from a fine-tuned BART model. |  Random Masking |  Custom Masking
---|---|---
Seed Text | does anyone have a travel discount code please? |  does anyone have a $<$mask$>$ code please?
Augmented Texts |  is anyone selling a travel discount logo please? do you have a travel discount app please? could anyone post a restaurant discount code please? |  does anyone have a discount code please? does anyone have a delivery code please? does anyone have a spare code please?

####  3.1.3 Quality and Throughput

For NLX, text quality depends on the developed grammar rules. A developer with
linguistic expertise can develop a precise grammar for the task. For highly
constrained domains, a grammar with acceptable coverage can be developed in a
matter of a few days. This is more challenging for complex and open domains,
such as LibriSpeech, and requires more efforts. Furthermore, while some
utterances tend to be syntactically correct, they may be less semantically
common. Below are examples of synthetic texts from NLX:

  * •

LibriSpeech: those hands could have been mexican

  * •

Photo Search: show me all pics of these this Tuesday

  * •

Q&A: search in burpee gardening for fences

For BART, as shown by Batra et. al. [16], random masking results in highly
diverse utterances. For custom masking, we get less variations as we mask only
nouns, verbs one at a time as compared to inserting random masks repeatedly N
times. Table 1 shows some samples generated from a fine-tuned BART model. The
BART method is free from the linguistic expertise and human intervention that
NLX does. It generates millions of high-quality utterances within 2 days using
a single GPU.

###  3.2 TTS

The TTS pipeline consists of a linguistic frontend and an acoustic backend.
The frontend takes plain text as input and emits its phonetic representation
(grapheme to phoneme) and an additional prosodic information for each token in
the output sequence. These are fed into the acoustic backend that consists of
a transformer-based prosody model which predicts phone-level F0 values and
duration, and a transformer-based spectral model that produces frame-level
mel-cepstral coefficients, F0 and periodicity features, and finally a sparse
WaveRNN-based neural vocoder that predicts the final waveform [20].

The acoustic model is trained on a multi-speaker dataset consisting of 170
hours of TTS-quality audio across 96 English speakers, while the WaveRNN
vocoder was trained for each speaker separately. We augment the 96 voices with
3 speeds and 5 pitch variations. This TTS system is fixed across all
experiments.

###  3.3 Audio Augmentation

The audio augmentation pipeline includes speed distortion to introduce speed
diversity by factors of 0.9, 1.0 and 1.1 for 40%, 20% and 40% of the data,
respectively. In addition, a pool of 1.72M noise tracks are used. We keep 40%
of audios noise-free, nearly 60% with one noise added, and $<$ 1% having 2-4
noises. The SNR has an N(12.50, 17.31) distribution.

###  3.4 ASR RNN-T Model

Our ASR model is a low-latency streamable RNN-T model, following the recipe in
[21]. Its audio encoder consists of 20 Emformer layers [22]. The RNN predictor
has 3 LSTM layers followed by a linear layer, and layer normalization in
between every two other layers. The joiner module is simply a ReLU activation
function followed by a single linear layer. We tokenize all texts to 5000
subword units [23] for LibriSpeech, and 4096 units for the two in-house
datasets. In total, our model has 78M trainable parameters for LibriSpeech,
and 76M for the in-house datasets. Adam optimizer is used for updating the
model parameters during training.

##  4 Experiments

We train ASR models for the three datasets described above. For each dataset,
we compare the baseline model trained with no data augmentation to models with
augmented TTS utterances. The TTS utterances are either created from the seed
text utterances, the rule-based NLX method, the neural BART method with either
random or custom masking, or a combination of the above. For all data
augmentation methods, the seed text is always included in the generated TTS
utterances.

###  4.1 Generated Data

The seed text data for the LibriSpeech, Photo Search and Q&A datasets is
comprised of 281k, 3350, and 100k utterances, respectively. Both NLX and BART
methods use an augmentation factors of 8.2x, 300x, and 8.0x, which resulted in
2.3M, 1.0M and 800k utterances respectively for the three datasets above.
Those were converted to speech utterances using the TTS system described in
Section 3.2. In order to simulate a situation of no human annotated utterances
from a target speech domain, all Photo Search and Q&A in-domain training data
is synthetically generated from our TTS system. When adding the synthetic
texts to the training corpus, we control the mixing ratio, i.e., how often
this data is seen by the model compared to other training data. We tune this
hyperparameter for each data and experiment separately. We found that higher
mixing ratios are acceptable for LibriSpeech, such that up to 50% of the data
seen is synthetic. For in-house datasets, we found that lower mixing ratios,
up to 10%, are optimal for best WER.

Table 2: Word error rate results (in percents) for baseline models and different data augmentation methods. | LibriSpeech | Photo Search | Q&A
---|---|---|---
Seed Data | 281 k | 3350 | 100 k
Augmented Data | 2.3 M | 1.0 M | 800 k
Test Set | test-other | test-clean | in-domain | in-domain-1 | in-domain-2
Baseline | 9.02 | 3.52 | 4.61 | 9.13 | 7.91
Seed Texts | 8.53 | 3.34 | 3.94 | 8.96 | 7.78
NLX | 8.65 | 3.31 | 4.18 | 8.34 | 7.15
BART C.M. | 8.51 | 3.24 | 3.94 | 8.41 | 7.15
BART R.M. | 8.14 | 3.13 | 3.94 | 8.32 | 7.07
NLX + BART R.M. | 8.36 | 3.14 | 4.03 | 8.08 | 6.78
NLX + BART R.M. + BART C.M. | 8.26 | 3.18 | 3.99 | 7.97 | 6.67

###  4.2 Results and Discussion

Table 2 lists the WERs for different trained models.

####  4.2.1 Seed Data without Augmentation

By introducing TTS from the seed texts only, the WERs for all the three
domains are reduced in different magnitudes, where Q&A sees a trivial relative
improvement ($<$ 2%), LibriSpeech achieves 5.2%-5.4%, and Photo Search
achieves 14.53%. Note that we only used 3350 utterances for the Photo Search
domain. Thus, for domain with simple and structured utterances, a few
thousands of text utterances are sufficient to achieve high performance ASR
with a WER $<$ 4%. TTS from seed texts is expected to improve both photo
search and Q&A, as the baseline model was trained without in-domain data.
Interestingly, we only observed the improvement on the simple photo search
domain.

####  4.2.2 Augmentation with BART or NLX

Adding synthetic texts to the seed data managed to improve WER results in many
but not all settings. Below we compare models trained with TA to models
trained with the seed texts only. NLX leads to a regression for both
LibriSpeech and the Photo Search domain, while a significant relative
improvement (6.1%-8.1%) for the Q&A domain is observed. For BART C.M., the WER
for the Photo Search domain has no regression, and sees improvements for the
other two domains. These improvements are further increased with BART R.M.

A few interesting insights arise from those comparisons. First, We see that
overall the BART neural approaches result in better WER gains compared to the
rule-based NLX. While this result may not be surprising, the superiority of
neural methods in TA followed by TTS for ASR data augmentation was not yet
verified in the existing literature. Second, we see that TA improves WER only
for the two datasets that have a large and diverse seed corpus, suggesting
that TA has a limited power when the seed data is small and structured.

####  4.2.3 Combining NLX and BART

We also experiment with the combination of synthetic data from different TA
approaches, where data is simply mixed together. We find that combining data
sources that result in gains individually (over just using the seed texts)
results in larger gains. For example, for the Q&A domain, all TA methods
improved over the seed data, and the ASR accuracy constantly improves from
BART R.M. to BART R.M. + NLX, to BART R.M. + NLX + BART C.M. However, for
LibriSpeech and the Photo Search domain, NLX does not result in large gains,
and we do not observe improvement after combining it with other data sources.

####  4.2.4 Effects of Augmented Factors

To further examine the effect of the size of augmented text corpus, Table 3
lists additional WERs from both test-clean and test-other where different
sizes of BART R.M. texts are introduced. From 2.3M to 7.3M, equivalent to an
augmentation factor of 8x to 26x, the WERs are fairly stable with a
fluctuation in relative improvement of less than 1%. We also augmented our
data by 103x and achieved WER improvements, where the test-clean set WER drops
from 3.11% to 3.00%. Similarly, for Q&A we also augmented the texts to 8.3M by
83x augmentation with BART R.M., but obtained no additional gains. Those
results indicate that in some cases TA methods may saturate and show no
further gains after a certain number of utterances are generated.

Note that synthetic data is generated at scale and at speed in this work.
Thus, we are less focused on small augmentation factors. Instead, for both Q&A
and photo search, we choose factors to introduce 4k hrs of audio or 2.7% of
the training data by duration. We also control TTS data used during training
by the mixing ratio. Both are to avoid overfitting to specific domains. For
LibriSpeech, we use a factor similar to the Q&A domain. The effect of TTS data
size on different domains requires further investigation in future work.

####  4.2.5 Other Factors

The size of the seed data plays a critical role, e.g. 10k (instead of 100k)
Q&A seed texts only achieves $<$ 5% relative improvement. Efforts to collect
high-quality seed text corpus are non-trivial, though dramatically reduced
comparing with collection of annotated paired audio-text data. Given the large
amounts of raw data available from external sources, e.g. web data, one might
argue that TA might not be needed. However, this type of data always requires
processing for privacy compliance, noise acceptance, and right distributions
for target domains. In this regard, it is highly desired to both increase
lexical diversities while reducing seed corpus size and to explore the
effective use of data from other resources.

In this work, audio augmentation with increased numbers of TTS voices and
audio distortions [4], e.g. noise addition and speed and pitch augmentation,
are beneficial for ASR improvement. Instead of 96 voices with audio distortion
used in this work, if we use only 12 voices and no audio distortion, the best
relative improvement that can be achieved is $<$ 4%. In future work, we will
be increasing voice diversity by enlarging our pool of TTS voices, and
enriching our customizable TTS with more personalized elements.

During the RNN-T decoding, external LMs are not used in this work, but only
the standard RNN-T beam search is employed. Nonetheless, comparing the
effectiveness of TTS vs. external LMs in utilizing the synthesized text is a
natural area for further exploration.

Table 3: LibriSpeech, Word Error Rate [%] vs. Text Size.

| 4.5 M | 5.9 M | 7.3 M | 29 M
---|---|---|---|---
Augmentation Factor | 16x | 21x | 26x | 103x
test-other | 8.21 | 8.18 | 8.21 | 8.10
test-clean | 3.10 | 3.11 | 3.11 | 3.00

##  5 Conclusion

In this work, we focus on text augmentation for speech data synthesis, which
we show to be an effective approach for improving ASR systems. Compared with
the rule-based method, a neural approach based on modern powerful pre-trained
language models is highly effective for synthetic text generation and the
subsequent ASR model training. Our results show that for the open-domain
LibriSpeech and Q&A datasets, neural text augmentation based on a fine tuned
BART model improves the relative performance by 9.76% - 15.68% compared to the
baselines. However, text augmentation results in no gain for simple domains
with highly-structured utterances. Furthermore, the additional WER gains by
simply mixing together data from different TA approaches confirmed synergy
among data from these techniques. This work suggests that it could be cost-
effective to improve an ASR system using an in-domain text corpus instead of
an expensive human-annotated audio corpus, and modern text generation
techniques should be considered as an important element in speech data
synthesis for improving ASR systems.

##  6 Acknowledgements

We would like to thank Xiao Yang, Jiale Zhi, Thilo Koehler, Didi Zhang, Azadeh
Nematzadeh, Qing He, Jilong Wu, and Tammy Stark for their help in data
collection, TTS, data augmentation, and helpful discussions and suggestions.

## References

  * [1] T. Hori, R. Astudillo, T. Hayashi, Y. Zhang, S. Watanabe, and J. L. Roux, ``Cycle-consistency training for end-to-end speech recognition,'' in _ICASSP_ , 2019.
  * [2] M. K. Baskar, L. Burget, S. Watanabe, R. F. Astudillo, and J. H. Černocký, ``Eat: Enhanced asr-tts for self-supervised speech recognition,'' in _ICASSP_ , 2021.
  * [3] R. Joshi and A. Singh, ``A simple baseline for domain adaptation in end to end asr systems using synthetic data,'' in _ECNLP_ , 2022, pp. 244–249.
  * [4] A. Fazel, W. Yang, Y. Liu, R. Barra-Chicote, Y. Meng, R. Maas, and J. Droppo, ``Synthasr: Unlocking synthetic data for speech recognition,'' in _Interspeech_ , 2021, pp. 896–900.
  * [5] C. Peyser, H. Zhang, T. N. Sainath, and Z. Wu, ``Improving performance of end-to-end ASR on numeric sequences,'' in _Interspeech_ , G. Kubin and Z. Kacic, Eds., 2019, pp. 2185–2189.
  * [6] R. Zevallos, N. Bel, G. Cámbara, M. Farrús, and J. Luque, ``Data augmentation for low-resource quechua asr improvement,'' _arXiv preprint arXiv:2207.06872_ , 2022.
  * [7] T.-Y. Hu, M. Armandpour, A. Shrivastava, J.-H. R. Chang, H. Koppula, and O. Tuzel, ``Synt++: Utilizing imperfect synthetic data to improve speech recognition,'' in _ICASSP_ , 2022, pp. 7682–7686.
  * [8] Z. Chen, A. Rosenberg, Y. Zhang, G. Wang, B. Ramabhadran, and P. J. Moreno, ``Improving speech recognition using gan-based speech synthesis and contrastive unspoken text selection.'' in _INTERSPEECH_ , 2020, pp. 556–560.
  * [9] A. Rosenberg, Y. Zhang, B. Ramabhadran, Y. Jia, P. Moreno, Y. Wu, and Z. Wu, ``Speech recognition with augmented synthesized speech,'' in _2019 IEEE automatic speech recognition and understanding workshop (ASRU)_ , 2019, pp. 996–1002.
  * [10] N. Rossenbach, M. Zeineldeen, B. Hilmes, R. Schlüter, and H. Ney, ``Comparing the benefit of synthetic training data for various automatic speech recognition architectures,'' in _2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_ , 2021, pp. 788–795.
  * [11] D. S. Park, W. Chan, Y. Zhang, C. Chiu, B. Zoph, E. D. Cubuk, and Q. V. Le, ``Specaugment: A simple data augmentation method for automatic speech recognition,'' in _Interspeech_ , G. Kubin and Z. Kacic, Eds., 2019, pp. 2613–2617.
  * [12] S. Ueno, M. Mimura, S. Sakai, and T. Kawahara, ``Data augmentation for asr using tts via a discrete representation,'' in _2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_ , 2021, pp. 68–75.
  * [13] B. Thai, R. Jimerson, D. Arcoraci, E. Prud'hommeaux, and R. Ptucha, ``Synthetic data augmentation for improving low-resource asr,'' in _2019 IEEE Western New York Image and Signal Processing Workshop (WNYISPW)_ , 2019, pp. 1–9.
  * [14] F. Biadsy, M. Ghodsi, and D. Caseiro, ``Effectively building tera scale maxent language models incorporating non-linguistic signals,'' in _Interspeech_ , 2017, pp. 2710–2714.
  * [15] M. Lewis, Y. Liu, N. Goyal, M. Ghazvininejad, A. Mohamed, O. Levy, V. Stoyanov, and L. Zettlemoyer, ``BART: Denoising sequence-to-sequence pre-training for natural language generation, translation, and comprehension,'' in _ACL_ , 2020, pp. 7871–7880.
  * [16] S. Batra, S. Jain, P. Heidari, A. Arun, C. Youngs, X. Li, P. Donmez, S. Mei, S. Kuo, V. Bhardwaj _et al._ , ``Building adaptive acceptability classifiers for neural nlg,'' in _EMNLP_ , 2021, pp. 682–697.
  * [17] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, ``Librispeech: an asr corpus based on public domain audio books,'' in _ICASSP_ , 2015, pp. 5206–5210.
  * [18] E. Reiter and R. Dale, _Building Natural Language Generation Systems_ , ser. Studies in Natural Language Processing. Cambridge University Press, 2000.
  * [19] S. Gururangan, A. Marasović, S. Swayamdipta, K. Lo, I. Beltagy, D. Downey, and N. A. Smith, ``Don't stop pretraining: Adapt language models to domains and tasks,'' in _ACL_ , 2020, pp. 8342–8360.
  * [20] C. Wu, Z. Xiu, Y. Shi, O. Kalinli, C. Fuegen, T. Koehler, and Q. He, ``Transformer-based acoustic modeling for streaming speech synthesis.'' in _Interspeech_ , 2021, pp. 146–150.
  * [21] D. Le, M. Jain, G. Keren, S. Kim, Y. Shi, J. Mahadeokar, J. Chan, Y. Shangguan, C. Fuegen, O. Kalinli, Y. Saraf, and M. L. Seltzer, ``Contextualized streaming end-to-end speech recognition with trie-based deep biasing and shallow fusion,'' in _Interspeech_ , 2021, pp. 1772–1776.
  * [22] Y. Shi, Y. Wang, C. Wu, C.-F. Yeh, J. Chan, F. Zhang, D. Le, and M. Seltzer, ``Emformer: Efficient memory transformer based acoustic model for low latency streaming speech recognition,'' in _ICASSP_ , 2021, pp. 6783–6787.
  * [23] T. Kudo and J. Richardson, ``Sentencepiece: A simple and language independent subword tokenizer and detokenizer for neural text processing,'' in _EMNLP_ , E. Blanco and W. Lu, Eds., 2018, pp. 66–71.


# Jovian Signal at BOREXINO

Saeed Ansarifard<EMAIL_ADDRESS>Yasaman Farzan<EMAIL_ADDRESS>School of physics, Institute for Research in Fundamental Sciences (IPM)
P.O.Box 19395-5531, Tehran, Iran

###### Abstract

The BOREXINO experiment has been collecting solar neutrino data since 2011,
providing the opportunity to study the variation of the event rate over a
decade. We find that at 96 % C.L., the rate of low energy events shows a time
modulation favoring a correlation with a flux from Jupiter. We present a new
physics model based on dark matter of mass 1-4 GeV captured by Jupiter that
can account for such modulation. We discuss how this scenario can be tested.

##  I Introduction

Accumulation of BOREXINO data on solar 7Be neutrino line over a decade [1]
provides the opportunity to study the time variation of the flux on Earth.
Obviously, because of the eccentricity of the Earth orbit,
$\epsilon=(1.6698\pm 0.0003)\%$
[2]111https://data.giss.nasa.gov/modelE/ar5plots/srorbpar.html the flux should
have a $\sim 3.3\%$ annual variation. Although the fluctuations in the
background of the experiment may contaminate this signal, thanks to the
installation of the temperature control and stabilization system, the
environmental conditions in the detector have been stabilized since 2019 [3].
There are also slowly decreasing radioactive backgrounds which have been
significantly reduced by the 2 final years of data taking [3].

Using the accurate BOREXINO data obtained after 2019, the annual variation of
7Be flux on Earth yields $\epsilon=(0.86\pm 0.62)\%$ which is about half the
actual eccentricity. On the other hand, with the first two years of data taken
during 2011$-$2013, the result tends to be 1.6 times the actual value:
$\epsilon=(2.74\pm 0.77)\%$. This tension cannot be attributed to the
fluctuations inside the Sun because the solar g-modes have much shorter
periodicity and average to zero over one year [4]. It is also worth noting
that the contribution of the cosmogenic ${}^{11}C$ due to the cosmic ray is
negligible in the energy range under study [1].

We entertain the possibility that the largest planet in the solar system,
Jupiter, may play a role. Indeed, we find that associating an event rate of
$\sim 1.5$ count per day per 100 tons of target $(\rm cpd/100t)$ to Jupiter,
the observed change in the modulation of events over years can be explained.
Such extra contribution can be accommodated within the average flux
uncertainty [5]. Measuring the direction of the flux over time, this
hypothesis can be tested. While dark matter with a mass of $1-4$ GeV can be
trapped inside Jupiter, it will evaporate from the Sun [6, 7, 8]. We build a
model based on this Jovian feature to explain an event rate of few per day at
BOREXINO associated to Jupiter. In our model, the annihilation of trapped DM
pairs produces a flux of milli-charged particles that leads to a signal at
BOREXINO mimicking the solar 7Be events.

##  II Hint for a signal from Jupiter

Figure 1: Relative modulation of the event rates at the detector due to the
orbital motion of the planets. The red (blue) curves shows the modulation
without (with) a contribution from Jupiter. The green area to the left
corresponds to the period when Jupiter is on the side of the perihelion of the
Earth orbit. The olive area to the right indicates the period for which
Jupiter is located closer to the aphelion. During this period, the background
conditions at the detector are most stable.

In the presence of background and a contribution from Jupiter, the
$\beta$-like event rate associated to the solar 7Be detection can be
parameterized as

$\mathcal{R}(t)=\frac{R_{sun}}{d^{2}_{sun}(t)}+\frac{R_{jup}}{d^{2}_{jup}(t)}+R_{B}$
(1)

where $d_{sun}$ ($d_{jup}$) is the distance from the Sun (Jupiter) to the
Earth in Astronomical Unit (AU) as a function of time $t$ which can be
obtained using Ref. [9]. $R_{sun}$ and $R_{jup}$ are respectively proportional
to the event rates induced by Sun and Jupiter. We take them to be constant in
time which is a reasonable assumption for the study of the annual variation.
Nevertheless, as discussed before $R_{B}$ at BOREXINO suffered a non-
negligible variation in the time period before 2019.

The temporal average of $\mathcal{R}(t)$ over years can be approximately
written as

$\langle\mathcal{R}\rangle\approx\dfrac{R_{sun}}{(1\ {\rm
AU})^{2}}+\dfrac{R_{jup}}{(5\ {\rm AU})^{2}}+\langle R_{B}\rangle$ (2)

As long as $R_{jup}\raise
1.29167pt\hbox{$\;<$\kern-7.5pt\raise-4.73611pt\hbox{$\sim\;$}}R_{sun}$, the
deviation of $\langle\mathcal{R}\rangle$ from the standard prediction will be
within uncertainties of the prediction [10]. Neglecting the time variation of
$R_{B}$, we define a quantity demonstrating the modulation in the event rate
due to varying $d_{sun}$ and $d_{jup}$

$\xi(t)\equiv\dfrac{\mathcal{R}(t)-\langle\mathcal{R}\rangle}{[R_{sun}/(1~{}{\rm
AU})^{2}]}.$ (3)

We expect when Jupiter is located closer to Earth aphelion (around 2020), the
annual variation of $\mathcal{R}(t)$ to be suppressed relative to the case
that Jupiter is on the opposite side (around 2013) because in the latter case
(when Jupiter is closer to perihelion), the annual maximums of
$R_{sun}/d_{sun}^{2}$ and $R_{jup}/d_{jup}^{2}$ add up. Fig. (1) shows the
behavior of $\xi(t)$ consistent with this expectation.

Utilizing the data from [1] we concentrate on the period Oct 2019 to Oct 2021
in which the temporal variation of $R_{B}$ is negligible [3]. Within standard
solar model and standard neutrino mass and mixing scheme, it is predicted
${R_{sun}}/{(1\ {\rm AU})^{2}}=25\pm 2\ (\rm cpd/100t)$ [11, 5]. Using Eq. (1)
and performing a Bayesian analysis [12, 13] with this prediction as an input,
we find the best fit value for $R_{jup}$ as

$\frac{R_{jup}}{(5\ {\rm AU})^{2}}=1.53^{+0.70}_{-0.81}\ (\rm cpd/100t)\ .$
(4)

The joint contour illustrated in Fig (2) shows an anti-correlation in the
$R_{B}-R_{jup}$ plane. As shown in the appendix, the likelihood ratio test
statistics disfavor the null Jovian signal (i.e., $R_{jup}=0$) at more than
$96\%\ \rm C.L$. Interestingly, the first two years of data taking which
corresponds to the expected excess of the annual modulation from Jupiter
favors this new solution by more than $2\sigma$. See the appendix for more
details.

Figure 2: $68\%$ and $95\%$ constraint contour on the free parameters
$R_{jup}$ and $R_{B}$. The analysis is performed using Eq. (1) and the last
two years of data from [1]. The predicted value of $R_{sun}$ is included with
a Gaussian prior.

##  III A model for Jovian Signal

In this section, we propose a toy model based on the absorption of Dark Matter
(DM) particles ($\chi$) in Jupiter and their subsequent annihilation to
intermediate particles that give rise to a signal at BOREXINO, mimicking that
of the solar 7Be neutrinos. We take the DM mass in the range 1-4 GeV [6, 7,
8]. As shown in [14], if $f_{\chi}=1\%$ of DM is composed of $\chi$, the spin-
dependent cross section of $\chi$ scattering off the protons can be as large
as $\sim 10^{-24}$ cm2. With a matter density of 1.3 ${\rm gr/cm^{3}}$ and
being 71 % composed of Hydrogen, the mean free path of $\chi$ on the Jovian
surface will be 2 cm. Thus, all $\chi$ particles reaching the Jupiter surface
will be absorbed, yielding an absorption rate of

$\mathcal{A}=4\pi R_{J}^{2}n_{\chi}v_{\chi}$

where $R_{J}=70000$ km is the Jupiter radius, $n_{\chi}=f_{\chi}(0.4~{}{\rm
GeV/cm^{3}})/m_{\chi}$ and $v_{\chi}\sim 10^{-3}c$ so

$\mathcal{A}=7\times 10^{25}~{}{\rm sec}^{-1}(f_{\chi}/0.01)({\rm
GeV}/m_{\chi}).$

The accumulated $\chi$ may annihilate to hadrons through the same coupling
responsible for scattering off nuclei. However, if we assume asymmetric
scenario where $n_{\chi}\gg n_{\bar{\chi}}$, this mode of annihilation will be
closed. Instead, we couple $\chi$ to a scalar $\phi$ with a mass $m_{\phi}\sim
20$MeV through the following scalar potential

$V=m_{\phi}^{2}|\phi|^{2}+(\lambda_{\phi\chi}|\phi|^{2}\chi^{2}+H.c.)\ .$ (5)

Notice that the global $U(1)$ symmetry under which $\chi\to\exp(-i\alpha)\chi$
is broken by the quadratic $\lambda_{\phi\chi}$ coupling but the $Z_{2}$
symmetry ($\chi\to-\chi$) remains unbroken and protects $\chi$ from decay. The
$\lambda_{\phi\chi}$ coupling can lead to $\chi\chi\to\phi\bar{\phi}$. For
$\lambda_{\phi\chi}<1$ and $m_{\chi}>1$ GeV

$\sigma_{\chi}=\sigma(\chi\chi\to\phi\bar{\phi})<3\times 10^{-29}\ {\rm
cm}^{2}.$

In order for $\chi\chi\to\phi\bar{\phi}$ to balance the $\chi$ capture in
Jupiter, the following relation should hold

$\mathcal{A}=2v_{\chi}^{J}(n_{\chi}^{J})^{2}V_{J}\sigma_{\chi}\ ,$

where $V_{J}=4\pi R_{J}^{3}/3$. $v_{\chi}^{J}$ and $n_{\chi}^{J}$ are
respectively the velocity and density of $\chi$ inside Jupiter. Taking the age
of Jupiter equal to $t_{J}=4.5\times 10^{9}$ years, we can write
$n_{\chi}^{J}\leq\mathcal{A}t_{J}/V_{J}$. In the limit that all captured
$\chi$ particles remain inside Jupiter the equality holds valid. The velocity
of $\chi$ should be between the thermal velocity in the Jupiter center ($\sim
10~{}{\rm km/sec}$) and the escape velocity (60 km$/$sec). Putting these
together, we find $\sigma_{\chi}\geq 10^{-36}$ cm2.

Annihilation $\phi\bar{\phi}\to\bar{\chi}\bar{\chi}$ at $T>$GeV can give rise
to nonzero $\bar{\chi}$ abundance but as long as we take $\sigma_{\chi}\gg
10^{-36}$ cm2, they annihilate away when $T\gg m_{\phi}$.

The $\phi$ pairs from annihilation of $\chi$ pairs will have an energy of
$m_{\chi}\sim$GeV. In order to lead to MeV scattering events at BOREXINO, the
$\phi$ particles have to first loose energy. Let us suppose there is a new
vector boson, $V$ with a mass of $(2-3)\times m_{\phi}$ and with couplings of
order of $g_{\phi}\sim O(1)$ to $\phi$ and of $g_{n}=O(10^{-3})$ to the
nucleons [15, 16]. Then,

$\sigma(\phi+N\to\phi+N)\sim\frac{g_{\phi}^{2}g_{n}^{2}}{4\pi m_{V}^{2}}\sim
10^{-32}{\rm cm}^{2}$

and

$\dfrac{\sigma(\phi+N\to\phi+N+V)}{\sigma(\phi+N\to\phi+N)}\sim\frac{g_{\phi}^{2}}{\pi^{2}}\left(\log\frac{\Delta
E}{m_{V}}\right)^{2}\sim{\rm few}\times 10^{-2}.$ (6)

With $g_{\phi}\sim 1$ and $g_{n}\sim 10^{-3}$, $\phi$ can go through a few
($\sim 5$) bremsstrahlung scatterings before leaving Jupiter. While in the two
body scattering, the average energy loss is $O(m_{V})$, in the bremsstrahlung
scattering the kinetic energy of $\phi$ will be shared between $V$ and final
$\phi$ and on average, the energy of $\phi$ can be reduced to half. $V$ will
immediately decay into a $\phi\bar{\phi}$ pair so a shower of $\phi$ and
$\bar{\phi}$ will be formed after $\sim 5$ bremsstrahlung scattering. The
energy of $\phi$ while leaving the planet will be therefore about
$2^{-5}m_{\chi}\sim m_{V}$. When the kinetic energy of $\phi$ drops below
$m_{V}$, the bremsstrahlung process stops. Moreover, in these energies
scattering off nuclei will be elastic and will not lead to a significant
energy loss. We therefore expect a flux of $\mathcal{A}(m_{\chi}/m_{V})$ of
$\phi$ particles with energies of $O(m_{V})$ streaming out Jupiter.

Figure 3: Schematic view of DM trapping inside Jupiter, its annihilation to
mediators and their cascading down to lower energies.

In principle, the $\phi$ particles can be milli-charged and interact with the
electrons inside the detector but then we expect a large signal induced by
$\phi$ coming from the Earth, too. To avoid such a contribution from the
Earth, we take $\phi$ to be neutral but assume that it decays to a pair of
milli-charged particles, $\phi\to C\bar{C}$ with a decay length much larger
than the earth radius. The $C$ particles with mass of $\sim 10$ MeV and charge
of $q\sim 10^{-4}e$ are still allowed by existing bounds [17, 18] and are even
favored by the EDGES data [19].

The scattering of $C$ off electrons inside BOREXINO with a recoil energy less
than $\sim 1$ MeV ($\Delta E\raise
1.29167pt\hbox{$\;<$\kern-7.5pt\raise-4.73611pt\hbox{$\sim\;$}}1$ MeV) can
mimic solar 7Be neutrino events. The cross section of this process can be
estimated as

$\sigma_{e}=\sigma(C+e\to C+e)|_{\Delta E<1~{}{\rm MeV}}\sim\frac{\alpha
q^{2}}{(m_{e}\Delta E)}$

The rate of electron scattering events per day per 100 tonnes with a recoil
energy less than $\sim$1 MeV is expected to be

$\mathcal{A}\times\frac{2m_{\chi}}{m_{V}}\times{\rm
day}\times\frac{\sigma_{e}}{4\pi d_{jup}^{2}}\times\frac{100~{}{\rm
t}}{2m_{p}}.$ (7)

The factor $2m_{\chi}/m_{V}$ takes care of the average multiplicity of the
$C+\bar{C}$ shower originating from a single $\chi$ annihilation. Thus, the
Jovian signal rate can be estimated as

$\frac{R_{jup}}{(5\ \rm AU)^{2}}\sim 1.5\frac{f_{\chi}}{3\times
10^{-4}}\left(\frac{q}{10^{-4}e}\right)^{2}\ \rm cpd/100t\ .$ (8)

##  IV Discussion and conclusions

We have analyzed the BOREXINO 7Be neutrino data collected from 2011 to 2021
and found that the time variation of the signal deviates from the expected
$1/d_{sun}^{2}$ behavior. In particular, the last time window of the data,
whose background was well under control, yields a value for the eccentricity
of the Earth orbit which is about half the known value. New physics scenarios
with 1 year periodicity such as pseudo-Dirac neutrinos [20, 21] are not able
to resolve this tension.

To explain the unexpected behavior, we have hypothesized a flux of new
particles from Jupiter that can induce a signal similar to that of the solar
7Be flux with a $1/d_{jup}^{2}$ time variation. We found that with such
contribution, the time variation of the signal at BOREXINO during 2019-2021
will be suppressed, explaining the lower extracted eccentricity. The best fit
for the average rate of events from Jupiter is about 6% of the 7Be events so
it can hide in the 6% uncertainty of the solar neutrino prediction. This
hypothesis can be tested by measuring the direction of incoming flux(es) over
time.

We have proposed a scenario that accounts for such signal from Jupiter, based
on three testable ingredients: (1) A fraction of $O(10^{-4})$ of dark matter
with mass of 1-4 GeV and with a strong spin-dependent scattering off nuclei
which can be tested by future dark matter search experiments [14]. Such dark
matter component can be trapped by Jupiter but not by Sun [6, 7, 8]. (2)
Intermediate particles, $\phi$ with mass of $O(10)$ MeV that are produced
inside Jupiter from dark matter annihilation and a light mediator ($V$) of
mass $O(10)$ MeV that couples to both $\phi$ and nuclei and leads to cascading
of the $\phi$ particles of energy of $\sim 1$ GeV inside Jupiter to a flux of
$\phi$ with kinetic energy of few MeV. The $V$ particles can be discovered by
the pion decay. (3) Milli-charged particles with mass of $O(10)$ MeV and
$q\sim 10^{-4}e$ (testable by DUNE [22]) that are produced by the decay of
$\phi$ on its way from Jupiter. The flux of milli-charged particle leads to a
signal at BOREXINO mimicking the solar 7Be line. Since the mediator is the
massless photon, the cross section increases at lower energies. Super-
Kamiokande further suppresses the Jovian signal by selecting only signals from
the Sun direction.

*

##  Appendix A Data Availability

All the material and the code for this article are publicly available in
Github

##  Appendix B Appendix: Hypothesis test statistics

In order to test whether the variation of data is because of the background
uncertainties, we randomly generate 30000 sample of simulated monthly binned
data with the Jovian signal, taking

$\\{\frac{R_{sun}}{(1\ \rm AU)^{2}}=25.3,\frac{R_{jup}}{(5\ \rm
AU)^{2}}=1.5,R_{B}=7.0\\}\ (\rm cpd/100t)$

and compare them with the standard case, taking

$\\{\frac{R_{sun}}{(1\ \rm AU)^{2}}=24.6,R_{jup}=0,R_{B}=9.2\\}\ (\rm
cpd/100t).$

Fig 4 shows the distribution of the $\boldsymbol{t}$ statistics which is
defined as

$\boldsymbol{t}=-2\log\left(\dfrac{\mathcal{L}({\rm SUN})}{\mathcal{L}({\rm
SUN+JUP})}\right)$ (9)

where $\mathcal{L}(F)$ is the likelihood of theory $F$ being true. The
BOREXINO result gives $\boldsymbol{t}=3.33$ which disfavors the standard
hypothesis with No-Jovian signal at $96.5\%$ C.L.

We repeat the test for the first two years of the data which yields the best
fit value

$\frac{R_{jup}}{(5\ \rm AU)^{2}}=1.8^{+0.8}_{-1.1}\rm(cpd/100t),$

consistent with the results from the last two years as shown in (4). In this
case because the variation of the background is not negligible, we have to use
the residual of the trend data provided in ref. [1]. The Jovian signal is
supported by more than $96\%$ C.L. Unsurprisingly, the rest of the BOREXINO
data taken during 2014 to 2019 is not able to discriminate the Jovian signal
due to the large uncertainties in the temporal variation of the background and
high correlation of the signal with the solar neutrino seasonal variation.

Figure 4: The distribution of the $\boldsymbol{t}$ statistics considering the
true signal comes only from Sun (red) or from both Sun and Jupiter (blue). The
dashed line shows the BOREXINO result.

##  Appendix C acknowledgments

This work has been supported by the European Union′s Framework Programme for
Research and Innovation Horizon 2020 under grant H2020-MSCA-
ITN-2019/860881-HIDDeN as well as under the Marie Sklodowska-Curie Staff
Exchange grant agreement No 101086085-ASYMMETRY. The authors would like to
acknowledge support from ICTP through the Associates Programme and from the
Simons Foundation through grant number 284558FY19. SA thanks A. M. Serenelli
for the helpful discussion. Fig. 2 is produced using GetDist package [23]. The
authors thank S. Jam for helping them to produce Fig. 3.

## References

  * Appel  _et al._ [2023] S. Appel _et al._ (BOREXINO), Astropart. Phys. 145, 102778 (2023), arXiv:2204.07029 [hep-ex] .
  * Simon  _et al._ [1994] J. L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze, G. Francou, and J. Laskar, AAP 282, 663 (1994).
  * Agostini  _et al._ [2020] M. Agostini _et al._ (BOREXINO), Nature 587, 577 (2020), arXiv:2006.15115 [hep-ex] .
  * Appourchaux  _et al._ [2010] T. Appourchaux _et al._ , Astron. Astrophys. Rev. 18, 197 (2010), arXiv:0910.0848 [astro-ph.SR] .
  * Vinyoles  _et al._ [2017] N. Vinyoles, A. M. Serenelli, F. L. Villante, S. Basu, J. Bergström, M. C. Gonzalez-Garcia, M. Maltoni, C. Peña Garay, and N. Song, Astrophys. J. 835, 202 (2017), arXiv:1611.09867 [astro-ph.SR] .
  * Li and Fan [2022] L. Li and J. Fan, JHEP 10, 186, arXiv:2207.13709 [hep-ph] .
  * French and Sher [2022] G. M. French and M. Sher, Phys. Rev. D 106, 115037 (2022), arXiv:2210.04761 [hep-ph] .
  * Buckley  _et al._ [2015] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, and G. Watt, Eur. Phys. J. C 75, 132 (2015), arXiv:1412.7420 [hep-ph] .
  * Rhodes [2019] B. Rhodes, Skyfield: High precision research-grade positions for planets and Earth satellites generator, Astrophysics Source Code Library, record ascl:1907.024 (2019), ascl:1907.024 .
  * Gonzalez-Garcia  _et al._ [2023] M. C. Gonzalez-Garcia, M. Maltoni, J. a. P. Pinheiro, and A. M. Serenelli, (2023), arXiv:2311.16226 [hep-ph] .
  * Esteban  _et al._ [2020] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09, 178, arXiv:2007.14792 [hep-ph] .
  * Torrado and Lewis [2019] J. Torrado and A. Lewis, Cobaya: Bayesian analysis in cosmology, Astrophysics Source Code Library, record ascl:1910.019 (2019), ascl:1910.019 .
  * Torrado and Lewis [2021] J. Torrado and A. Lewis, JCAP 05, 057, arXiv:2005.05290 [astro-ph.IM] .
  * Li  _et al._ [2023] Y. Li, Z. Liu, and Y. Xue, JCAP 05, 060, arXiv:2209.04387 [hep-ph] .
  * Gninenko and Krasnikov [1998] S. N. Gninenko and N. V. Krasnikov, Phys. Lett. B 427, 307 (1998), arXiv:hep-ph/9802375 .
  * Altegoer  _et al._ [1998] J. Altegoer _et al._ (NOMAD), Phys. Lett. B 428, 197 (1998), arXiv:hep-ex/9804003 .
  * Prinz  _et al._ [1998] A. A. Prinz _et al._ , Phys. Rev. Lett. 81, 1175 (1998), arXiv:hep-ex/9804008 .
  * Arefyeva  _et al._ [2022] N. Arefyeva, S. Gninenko, D. Gorbunov, and D. Kirpichnikov, Phys. Rev. D 106, 035029 (2022), arXiv:2204.03984 [hep-ph] .
  * Berlin  _et al._ [2019] A. Berlin, N. Blinov, G. Krnjaic, P. Schuster, and N. Toro, Phys. Rev. D 99, 075001 (2019), arXiv:1807.01730 [hep-ph] .
  * Ansarifard and Farzan [2023] S. Ansarifard and Y. Farzan, Phys. Rev. D 107, 075029 (2023), arXiv:2211.09105 [hep-ph] .
  * Anamiati  _et al._ [2018] G. Anamiati, R. M. Fonseca, and M. Hirsch, Phys. Rev. D 97, 095008 (2018), arXiv:1710.06249 [hep-ph] .
  * Magill  _et al._ [2019] G. Magill, R. Plestid, M. Pospelov, and Y.-D. Tsai, Phys. Rev. Lett. 122, 071801 (2019), arXiv:1806.03310 [hep-ph] .
  * Lewis [2019] A. Lewis,  (2019), arXiv:1910.13970 [astro-ph.IM] .