# On the Unintended Social Bias of Training
Language Generation Models with Data from
Local Media

Omar U. Florez
Conversational AI Research
Capital One
<EMAIL_ADDRESS>

###### Abstract

There are concerns that neural language models may preserve some of the
stereotypes of the underlying societies that generate the large corpora needed
to train these models. For example, gender bias is a significant problem when
generating text, and its unintended memorization could impact the user
experience of many applications (e.g., the smart-compose feature in Gmail).

In this paper, we introduce a novel architecture that decouples the
representation learning of a neural model from its memory management role.
This architecture allows us to update a memory module with an equal ratio
across gender types addressing biased correlations directly in the latent
space. We experimentally show that our approach can mitigate the gender bias
amplification in the automatic generation of articles news while providing
similar perplexity values when extending the Sequence2Sequence architecture.

##  1 Introduction

Neural Networks have proven to be useful for automating tasks such as question
answering, system response, and language generation considering large textual
datasets. In learning systems, bias can be defined as the negative
consequences derived by the implicit association of patterns that occur in a
high-dimensional space. In dialogue systems, these patterns represent
associations between word embeddings that can be measured by a Cosine distance
to observe male- and female-related analogies that resemble the gender
stereotypes of the real world. We propose an automatic technique to mitigate
bias in language generation models based on the use of an external memory in
which word embeddings are associated to gender information, and they can be
sparsely updated based on content-based lookup.

The main contributions of our work are the following:

  * •

We introduce a novel architecture that considers the notion of a Fair Region
to update a subset of the trainable parameters of a Memory Network.

  * •

We experimentally show that this architecture leads to mitigate gender bias
amplification in the automatic generation of text when extending the
Sequence2Sequence model.

##  2 Memory Networks and Fair Region

As illustrated in Figure 1, the memory $M$ consists of arrays $K$ and $V$ that
store addressable keys (latent representations of the input) and values (class
labels), respectively as in Kaiser et al. (2017). To support our technique, we
extend this definition with an array $G$ that stores the gender associated to
each word, e.g., actor is male, actress is female, and scientist is no-gender.
The final form of the memory module is as follows:

$M=(K,V,G).$

Figure 1: A Fair Region in memory $M$ consists of the most similar keys to $h$
given a uniform distribution over genders (e.g., $1$: male, and $2$: female).
The input consists of a sequence tokens annotated with gender information,
e.g., (The, 0), (president, 0), (gave, 0), (her, 2), (speech, 0).

A neural encoder with trainable parameters $\theta$ receives an observation
$x$ and generates activations $h$ in a hidden layer. We want to store a
normalized $h$ (i.e., $\left\lVert h\right\rVert=1$) in the long-term memory
module $M$ to increase the capacity of the encode. Hence, let $i_{max}$ be the
index of the most similar key

$i_{max}=argmax_{i}\\{h\cdot K[i]\\},$

then writing the triplet $(x,y,g)$ to $M$ consist of:

$\displaystyle K[i_{max}]$ $\displaystyle=\left\lVert
h+K[i_{max}]\right\rVert$ $\displaystyle V[i_{max}]$ $\displaystyle=y$
$\displaystyle G[i_{max}]$ $\displaystyle=g.$

However, the number of word embeddings does not provide an equal
representation across gender types because context-sensitive embeddings are
severely biased in natural language, Zhao et al. (2017). For example, it has
been shown in that man is closer to programmer than woman, Bolukbasi et al.
(2016). Similar problems have been recently observed in popular work embedding
algorithms such as Word2Vec, Glove, and BERT, Kurita (2019).

We propose the update of a memory network within a Fair Region in which we can
control the number of keys associated to each particular gender. We define
this region as follows.

######  Definition 2.1.

(Fair Region) Let $h$ be an latent representation of the input and $M$ be an
external memory. The male-neighborhood of $h$ is represented by the indices of
the $n$-nearest keys to $h$ in decreasing order that share the same gender
type male as $\\{i^{m}_{1},...,i^{m}_{k}\\}=KNN(h,n,male)$. Running this
process for each gender type estimates the indices $i^{m}$, $i^{f}$, and
$i^{ng}$ which correspond to the male, female, and non-gender neighborhoods.
Then, the FairRegion of $M$ given $h$ consists of $K[i^{m};i^{f};i^{ng}]$.

The Fair Region of a memory network consists of a subset of the memory keys
which are responsible for computing error signals and generating gradients
that will flow through the entire architecture with backpropagation. We do not
want to attend over all the memory entries but explicitly induce a uniform
gender distribution within this region. The result is a training process in
which gender-related embeddings equally contribute in number to the update of
the entire architecture. This embedding-level constraint prevents the
unconstrained learning of correlations between a latent vector $h$ and similar
memory entries in $M$ directly in the latent space considering explicit gender
indicators.

##  3 Language Model Generation

Figure 2: Architecture of the neural dialogue model that incorporates the
memory approach discussed in Section 1. The figure shows the $i^{th}$ decoding
step of the word $\hat{y_{i}}$ given the sparse update within a Fair Region
centered at $h^{deco}$.

Our goal is to leverage the addressable keys of a memory augmented neural
network and the notion of fair regions discussed in Section2 to guide the
automatic generation of text. Given an encoder-decoder architecture Sutskever
et al. (2014); Bahdanau et al. (2015), the inputs are two sentences $x$ and
$y$ from the source and target domain, respectively. An LSTM encoder outputs
the context-sensitive hidden representation $h^{enco}$ based on the history of
sentences and an LSTM decoder receives both $h^{enco}$ and $y$ and predicts
the sequence of words $\hat{y}$. At every timestep of decoding, the decoder
predicts the $i^{th}$ token of the output $\hat{y}$ by computing its
corresponding hidden state $h^{deco}_{i}$ applying the recurrence

$\displaystyle h^{deco}_{i}$ $\displaystyle=LSTM(y_{i-1},h^{deco}_{i-1}).$

Instead of using the decoder output $h_{i}^{deco}$ to directly predict the
next word as a prediction over the vocabulary $O$, as in Miller et al. (2016).
We combine this vector with a query to the memory module to compute the
embedding vector $h^{fair}_{i}$. We do this by computing an attention score
Bahdanau et al. (2015) with each key of a Fair Region. The attention logits
become the unormalized probabilities of including their associated values for
predicting the $i^{th}$ token of the response $\hat{y}$. We then argmax the
most likely entry in the output vocabulary $O$ to obtain the $i^{th}$
predicted token of the response $\hat{y}$. More formally,

$\displaystyle K^{fair}$ $\displaystyle=FairRegion(h^{deco},K,n)$
$\displaystyle\alpha_{i}$ $\displaystyle=Softmax(h^{deco}_{i}\cdot K^{fair})$
$\displaystyle h_{i}^{fair}$ $\displaystyle=Wtanh(\alpha_{i}\cdot K^{fair})$
$\displaystyle\hat{y}$ $\displaystyle=Softmax(h_{i}^{fair})$
$\displaystyle\hat{y}_{i}$ $\displaystyle=O[argmax_{j}(\hat{y}[j])].$

Naturally, the objective function is to minimize the cross entropy of actual
and generated content:

$J(\theta)=-\sum_{j=1}^{N}\sum_{i=1}^{m}y_{i}^{j}\ log\ p(\hat{y}_{i}^{j})$

where $N$ is the number of training documents, $m$ indicates the number of
words in the generated output, and $y_{i}^{j}$ is the one-hot representation
of the $i^{th}$ word in the target sequence.

##  4 Bias Amplification

As originally introduced by Zhao et al. (2017), we compute the bias score of a
word $x$ considering its word embedding $h^{fair}(x)$111For Seq2Seq neural
models, this word embedding is the output of the decoder component
$h^{deco}(x)$ and two gender indicators (words man and woman). For example,
the bias score of scientist is:

$b(scientist,man)=\frac{\left\lVert
h^{fair}(scientist)\right\rVert\cdot\left\lVert
h^{fair}(man)\right\rVert}{\left\lVert h^{fair}(scientist)\cdot
h^{fair}(man)+h^{fair}(scientist)\cdot h^{fair}(woman)\right\rVert}.$

If the bias score during testing is greater than the one during training,

$b^{test}(scientist,man)-b^{train}(scientist,man)>0,$

then the bias of man towards scientist has been amplified by the model while
learning such representation, given training and testing datasets similarly
distributed.

| Perplexity | Bias Amplification
---|---|---
Model | All | Peru | Mexico | All | Peru-Mexico | Mexico-Peru
Seq2Seq | 13.27 | 15.31 | 15.61 | +0.18 | +0.25 | +0.21
Seq2Seq+Attention | 10.73 | 13.25 | 14.08 | +0.25 | +0.32 | +0.29
SeqSeq+FairRegion | 10.79 | 13.04 | 13.91 | +0.09 | +0.17 | +0.15
Table 1: Perplexity and Bias Amplification results on the datasets of crawled
newspapers.

##  5 Experiments

###  5.1 Dataset

We evaluate our proposed method in datasets crawled from the websites of three
newspapers from Chile, Peru, and Mexico.

To enable a fair comparison, we limit the number of articles for each dataset
to 20,000 and the size of the vocabulary to the 18,000 most common words.
Datasets are split into 60%, 20%, and 20% for training, validation, and
testing. We want to see if there are correlations showing stereotypes across
different nations. Does the biased correlations learned by an encoder transfer
to the decoder considering word sequences from different countries?

###  5.2 Baselines

We compare our approach Seq2Seq+FairRegion, an encoder-decoder architecture
augmented with a Fair Region, with the following baseline models:

  1. 1.

Seq2Seq Sutskever et al. (2014): An encoder-decoder architecture that maps
between sequences with minimal assumptions on the sequence structure and that
is able to remember long term dependencies by mapping the source sentence into
a fixed-length vector.

  2. 2.

Seq2Seq+Attention Bahdanau et al. (2015): Similar to Seq2Seq, this
architecture automatically attends to parts of the input that can be relevant
to predict the target word.

###  5.3 Training Settings

For all the experiments, the size of the word embeddings is 256. The encoders
and decoders are bidirectional LSTMs of 2-layers with state size of 256 for
each direction. For the Seq2Seq+FairRegion model, the number of memory entries
is 1,000. We train all models with Adam optimizer Kingma & Ba (2014) with a
learning rate of $0.001$ and initialized all weights from a uniform
distribution in $[-0.01,0.01]$. We also applied dropout Srivastava et al.
(2014) with keep probability of $95.0\%$ for the inputs and outputs of
recurrent neural networks.

###  5.4 Fair Region Results in Similar Perplexity

We evaluate all the models with test perplexity, which is the exponential of
the loss. We report in Table 1 the average perplexity of the aggregated
dataset from Peru, Mexico, and Chile, and also from specific countries.

Our main finding is that our approach (Seq2Seq+FairRegion) shows similar
perplexity values ($10.79$) than the Seq2Seq+Attention baseline model
($10.73$) when generating word sequences despite using the Fair Region
strategy. These results encourage the use of a controlled region as an
automatic technique that maintains the efficacy of generating text. We observe
a larger perplexity for country-based datasets, likely because of their
smaller training datasets.

###  5.5 Fair Region Controls Bias Amplification

We compute the bias amplification metric for all models, as defined in Section
4, to study the effect of amplifying potential bias in text for different
language generation models.

Table 1 shows that using Fair Regions is the most effective method to mitigate
bias amplification when combining all the datasets (+0.09). Instead, both
Seq2Seq (+0.18) and Seq2Seq+Attention (+0.25) amplify gender bias for the same
corpus. Interestingly, feeding the encoders with news articles from different
countries decreases the advantage of using a Fair Region and also amplifies
more bias across all the models. In fact, training the encoder with news from
Peru has, in general, a larger bias amplification than training it with news
from Mexico. This could have many implications and be a product of the writing
style or transferred social bias across different countries. We take its
world-wide study as future work.

##  6 Conclusions

Gender bias is an important problem when generating text. Not only smart
composer or auto-complete solutions can be impacted by the encoder-decoder
architecture, but the unintended harm made by these algorithms could impact
the user experience in many applications. We also show the notion of bias
amplification applied to this dataset and results on how bias can be
transferred between country-specific datasets in the encoder-decoder
architecture.

## References

  * Bahdanau et al. (2015) Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio.  Neural machine translation by jointly learning to align and translate.  In _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.  URL http://arxiv.org/abs/1409.0473.
  * Bolukbasi et al. (2016) Tolga Bolukbasi, Kai-Wei Chang, James Zou, Venkatesh Saligrama, and Adam Kalai.  Man is to computer programmer as woman is to homemaker? debiasing word embeddings.  NIPS’16, pp. 4356–4364, 2016.
  * Kaiser et al. (2017) Lukasz Kaiser, Ofir Nachum, Aurko Roy, and Samy Bengio.  Learning to remember rare events.  2017\.  URL https://openreview.net/pdf?id=SJTQLdqlg.
  * Kingma & Ba (2014) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  CoRR, 2014.
  * Kurita (2019) Keita Kurita.  Measuring bias in contextualized word representations.  In _1st ACL Workshop on Gender Bias for Natural Language Processing 2019_ , pp. 166–172, 01 2019.  doi: 10.18653/v1/W19-3823.
  * Miller et al. (2016) Alexander Miller, Adam Fisch, Jesse Dodge, Amir-Hossein Karimi, Antoine Bordes, and Jason Weston.  Key-value memory networks for directly reading documents.  In _Proceedings of the 2016 Conference on Empirical Methods in Natural Language Processing_ , pp. 1400–1409, Austin, Texas, November 2016. Association for Computational Linguistics.  doi: 10.18653/v1/D16-1147.  URL https://www.aclweb.org/anthology/D16-1147.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _J. Mach. Learn. Res._ , 15(1):1929–1958, January 2014.  ISSN 1532-4435.  URL http://dl.acm.org/citation.cfm?id=2627435.2670313.
  * Sutskever et al. (2014) Ilya Sutskever, Oriol Vinyals, and Quoc V Le.  Sequence to sequence learning with neural networks.  In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger (eds.), _Advances in Neural Information Processing Systems 27_ , pp. 3104–3112. Curran Associates, Inc., 2014.  URL http://papers.nips.cc/paper/5346-sequence-to-sequence-learning-with-neural-networks.pdf.
  * Zhao et al. (2017) Jieyu Zhao, Tianlu Wang, Mark Yatskar, Vicente Ordonez, and Kai-Wei Chang.  Men also like shopping: Reducing gender bias amplification using corpus-level constraints.  In _Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , 2017.


11institutetext: JD Technology Group
11email: {songwei11,yueyanghao,zhangyajie23,zhangzhengchen1,
<EMAIL_ADDRESS>
# Multi-Speaker Multi-Style Speech Synthesis with Timbre and Style
Disentanglement††thanks: This work was supported by the National Key R&D
Program of China under Grant No. 2020AAA0108600.

Wei Song  Yanghao Yue  Ya-jie Zhang  Zhengchen Zhang  Youzheng Wu  Xiaodong He

###### Abstract

Disentanglement of a speaker’s timbre and style is very important for style
transfer in multi-speaker multi-style text-to-speech (TTS) scenarios. With the
disentanglement of timbres and styles, TTS systems could synthesize expressive
speech for a given speaker with any style which has been seen in the training
corpus. However, there are still some shortcomings with the current research
on timbre and style disentanglement. The current method either requires
single-speaker multi-style recordings, which are difficult and expensive to
collect, or uses a complex network and complicated training method, which is
difficult to reproduce and control the style transfer behavior. To improve the
disentanglement effectiveness of timbres and styles, and to remove the
reliance on single-speaker multi-style corpus, a simple but effective timbre
and style disentanglement method is proposed in this paper. The FastSpeech2
network is employed as the backbone network, with explicit duration, pitch,
and energy trajectory to represent the style. Each speaker’s data is
considered as a separate and isolated style, then a speaker embedding and a
style embedding are added to the FastSpeech2 network to learn disentangled
representations. Utterance level pitch and energy normalization are utilized
to improve the decoupling effect. Experimental results demonstrate that the
proposed model could synthesize speech with any style seen during training
with high style similarity while maintaining very high speaker similarity.

###### Keywords:

Speech Synthesis Style Transfer Disentanglement.

##  1 Introduction

With the development of deep learning technology in the last decade, speech
synthesis technology has evolved from traditional statistics-based speech
synthesis [1] to end-to-end based [2, 3, 4, 5, 6] and made great advancements.
The current speech synthesis technology has been able to synthesize speech
with high naturalness and high fidelity, and even some research [7] has been
able to synthesize speech that human beings cannot distinguish between true
recordings.

Although the great achievement in speech synthesis, there still exists a large
improvement room for expressive speech synthesis, especially for multi-speaker
multi-style Text to Speech (TTS) with cross-speaker style transfer [8, 9, 10].
Synthesizing speech with a target speaker’s timbre and other speaker’s style
could further increase the application scenarios and expressiveness of the TTS
system.

In order to synthesize more expressive speech, some researches [11, 8, 10] do
style decoupling by using a single-speaker multi-style recording data to learn
the style representation. However, it is difficult and expensive to collect
these kinds of data in lots of scenarios. Some other works [12, 9] try to
remove this data restriction by learning a decoupled style representation, but
most of them use complex network structure or rely on complicated training
method, which makes it difficult to re-produce the experimental results or
infeasible to control the network performance during inference.

A jointly trained reference encoder is used to learn implicit style
representation in [12, 13]. After the model is trained, audio with a different
style or even from a different speaker could be taken as the reference audio
to synthesize speech with the desired style while keeping the timbre
unchanged. However, the reference-based method is unstable which usually
generates unexpected style, and it’s non-trivial to choose the reference
audio.

Pan et al. [8] use prosody bottle-neck feature to learn a compact style
representation, and paper [10] uses a separate style encoder and speaker
encoder to disentangle the speaker’s timbre and style, and cycle consistent
loss is used to improve the disentanglement effect, a complex neural network
and complicated training objectives are required to achieve good performance.
Both [10] and [8] require single-speaker multi-style corpus to learning the
disentanglement, which restricts the flexibility of their proposed models.

In this paper, we propose a simple but effective expressive speech synthesis
network that disentangles speakers’ timbres and styles, which makes it
available to do multi-speaker multi-style speech synthesis. The proposed
system does not require a multi-speaker multi-style corpus, each speaker’s
data is considered as an isolated style. The proposed network is similar to
work [9], but we use FastSpeech2 [4] as the network backbone, which removes
the skip/repeat pronunciation issue caused by the attention mechanism. Prosody
features (duration, pitch, and energy) are used directly to improve the
disentanglement effectiveness of timbre and style, furthermore, utterance
level pitch and energy normalization (UttNorm) are used to prevent identity
leakage from prosody features.

##  2 The Proposed Model

The proposed network utilizes FastSpeech2 [4] as the network backbone and
applies utterance level pitch and energy normalization to achieve a better
decoupling effect. MelGAN [14] is used as the neural vocoder to convert
acoustic features to speech.

###  2.1 The Network Structure Of Proposed Network

FastSpeech2 [4] network architecture is used as the network backbone for the
proposed model, which consists of a phoneme encoder to learn syntactic and
semantic features, a mel-spectrogram [3] decoder to generate frame-level
acoustic features, and a variance adaptor to learn style-related features. The
network structure is shown in Fig. 1 (a).

A style embedding is introduced to the network to learn a style-dependent
variance adaptor, the network structure of variance adapter is illustrated in
Fig. 1 (b) and Fig. 1 (c). To decouple timbre and style, the speaker embedding
is moved to the input of the decoder to make the decoder timbre dependent.

One important purpose of the proposed model is to remove the dependency of the
single-speaker multi-style corpus, so each speaker’s corpus is considered as a
unique style and we could learn style representation from other speakers’
corpus. Actually, during the training procedure, the speaker id and style id
is identical, then it is very important to prevent style embedding from
leaking into the backbone network, which means style embedding should only be
used for style feature prediction and never be exposed to the backbone
network. So in this proposed network, style embedding is only used in the
variance adaptor to ensure that the whole network learns the speaker’s timbre
by speaker embedding instead of style embedding, and style is only affected by
style embedding.

Figure 1: An illustration of the proposed network architecture. (a) Network
structure of the acoustic model. (b) The variance adaptor structure with style
embedding as conditional input. (c) The structure of variance predictor.

###  2.2 Utterance Level Feature Normalization

The speaking style is represented in many aspects, such as the duration of
each syllable, the fundamental frequency (F0), and the trajectory of pitch and
energy. However, these scalars could still contain speaker identity and impact
the decoupling effect. For example, female speakers usually have higher F0
[15] than male speakers, and excited speakers express higher energy value, so
the pitch and energy features contain speaker timbre information to some
extent, then it’s essential to normalize the style features to disentangle a
speaker’s timbre and style.

In this paper, instead of speaker level normalization (SpkNorm), utterance
level pitch and energy normalization (UttNorm) are used to remove speaker
identity from the style features for better timbre and style disentanglement.
The style difference in each utterance will lead to a statistics difference
between speaker level and utterance level statistics, this difference could
cause identity leakage and affect the decoupling effect, especially for the
unprofessional, noisy recordings. UttNorm could eliminate this statistic
difference and improve the timbre and style disentanglement effectiveness.

##  3 Experimental Setup

Three open-sourced Chinese mandarin corpora and three internal Chinese
mandarin corpora with distinctive styles are used to train the proposed model.
The open-sourced corpus includes CSMSC 111https://www.data-
baker.com/open˙source.html, which is recorded by a female speaker, and the
MST-Originbeat
222http://challenge.ai.iqiyi.com/detail?raceId=5fb2688224954e0b48431fe0 [16],
which consists of recordings from a female speaker and a male speaker. Details
of the training data are listed in table 1.

Table 1: The description of the dataset used in the experiment. There are three open-sourced corpus (with a * in the speaker name) and three internal corpora. Speaker | #Utterance | Style | Gender
---|---|---|---
CSMSC* | 10000 | normal |  female
Originbeat-S1* | 5000 | normal | female
Originbeat-S2* | 5000 | normal | male
C1 | 10000 | children story | female
F1 | 450 |  news-broadcasting | female
M1 | 35000 | story-telling | male

The training data waves are converted to 16kHz, 16bit depth, and then scaled
to 6dB in our experiments. The extracted phoneme labels and processed speech
waves are aligned by the MFA [17] tool to detect the phoneme boundary. The
80-band mel-scale spectrogram is extracted as the training target with a 12ms
hop size and 48ms window size. Pitch is extracted by using the PyWORLD
333https://github.com/JeremyCCHsu/Python-Wrapper-for-World-Vocoder toolkit.
Both pitch and energy trajectories are normalized in utterance level to remove
speaker identity.

The proposed model is trained by Adam [18] optimizer with a batch size of 32
and noam [19] learning rate schedule. The learning rate is warmed up to a
maximum value of 1e-3 in the first 4000 steps and then exponentially decayed.
The model is trained by 400,000 steps and the network is regularized by weight
decay with a weight of 1e-6.

##  4 Experimental Results

This chapter shows the experimental results of the proposed timbre and style
disentanglement network. We encourage the readers to listen to the synthesized
speeches on our demo page
444https://weixsong.github.io/demos/MultiSpeakerMultiStyle/.

###  4.1 Subjective Evaluation

To evaluate the proposed method, Mean Opinion Score (MOS) evaluation is
conducted to evaluate the speaker similarity and style similarity for the
synthesized speech of different speakers and different styles. The target
speakers are from the 3 open-sourced corpora, and our internal news-
broadcasting (F1), children story (C1), and story-telling (M1) styles are used
as the target style. Twenty utterances are synthesized for each speaker and
style combination, listened to by 15 testers. During each MOS evaluation, the
speaker similarity or style similarity is the only point that testers need to
focus on.

Table 2: The 5 points MOS results with a confidence interval of 95%. 5 means the style or the timbre is exactly the same as the reference, and 1 means totally different. Speaker | Speaker Similarity | Style Similarity
---|---|---
C1 | F1 | M1 | C1 | F1 | M1
CSMSC | 4.50$\pm$0.05 | 4.47$\pm$0.06 | 4.49$\pm$0.04 | 4.31$\pm$0.03 | 4.62$\pm$0.03 | 4.11$\pm$0.06
Originbeat-S1 | 4.36$\pm$0.08 | 4.66$\pm$0.02 | 4.29$\pm$0.06 | 4.13$\pm$0.03 | 4.18$\pm$0.05 | 4.08$\pm$0.08
Originbeat-S2 | 4.36$\pm$0.07 | 4.59$\pm$0.04 | 4.23$\pm$0.07 | 4.23$\pm$0.04 | 4.30$\pm$0.03 | 4.10$\pm$0.08

MOS results of speaker similarity and style similarity are shown in Table 2,
we could see that for a given target speaker, the speaker similarity of
speeches in different style are very high and consistent, and the style
similarity is also very high with the target style, which indicates that the
proposed model achieves excellent disentanglement effect for speaker’s timbre
and style.

###  4.2 Ablation Study Of Utterance Level Feature Normalization

To verify the effectiveness of UttNorm, another female corpus(112 sentences)
crawled from a podcast, even with some background noise is chosen to
demonstrate the impact.

In our experiments, due to the large variety of this noisy corpus from the
oral podcast, a speaker similarity MOS evaluation and a preference evaluation
are conducted to evaluate the performance of different normalization methods.
We use speaker F1 as the target speaker and the style from this noisy corpus
as the target style, then 20 sentences in the podcast domain are synthesized
and listened to by 15 testers.

According to the MOS results for speaker similarity, when UttNorm is used, the
speaker similarity increases from 3.82$\pm$0.05(SpkNorm) to
3.91$\pm$0.06(UttNorm), which showed that UttNorm could facilitate the
decoupling effect of timbre and style for data with high variance, especially
found data, noisy data or podcast data.

From the AB preference evaluation results, the proportions of preference are
0.17(SpkNorm), 0.35(No Preference), 0.48(UttNorm), and the p-value is less
than 0.001. It could be found that the listeners strongly prefer the results
from UttNorm, which proves the effectiveness of utterance normalization.

###  4.3 Demonstration Of The Proposed Model

Synthesized speeches with a given speaker and different styles are shown in
Fig. 2 (a), and speeches for a given style and different speakers are shown in
Fig. 2 (b). The fundamental frequency is consistent for the three different
styles in Fig. 2 (a), and all the speakers’ pitch trajectories follow
basically the same curve but with different fundamental frequency in Fig. 2
(b), this explains the decoupling effectiveness to some extent.

Figure 2: A mel-spectrogram illustration of the proposed multi-speaker multi-
style TTS model. (a) Speaker F1 is selected as the target speaker, the target
style is the original news-broadcasting style, children story style, and
story-telling style. (b) News-broadcasting style is selected as target style,
the speakers are CSMSC, Originbeat-S1, and Originbeat-S2 respectively, from
top to bottom.

###  4.4 Style Transition Illustration

As both a speaker ID and a style ID should be sent to the network to
synthesize the target speaker’s speech with a given style, we could also use
the speaker ID to generate a style embedding that represents the source style,
and the embedding from the style ID represents the target style. By combining
the style embeddings from source and target with different weights, we could
generate speech with different target style intensities.

To demonstrate the continuity of the learned style embedding representation, a
style transition example is given in Fig. 3. From the pitch trajectory in each
synthesized speech we could find that the fundamental frequency is identical
for different target style weight, indicating that the proposed model keep the
timbre unchanged when synthesizing different style speech and achieves
outstanding timbre and style disentanglement effect.

Figure 3: A style transition illustration. Speaker F1 is selected as the
target speaker and style M1 is selected as the target style, this figure shows
the gradual style transition from style F1 to style M1 from top to bottom,
with different target style weights.

##  5 Conclusions

A simple but effective speaker’s timbre and style disentanglement network is
proposed in this paper, which eliminates the reliance on a single-speaker
multi-style corpus. The proposed network learns a style-dependent variance
adaptor and a speaker-dependent mel-spectrogram prediction decoder. The style-
related features are predicted by the variance adaptor with the guidance of
style embedding, while the timbre is learned by the mel-spectrogram decoder
with the control of speaker embedding. Utterance level feature normalization
is proposed to prevent speaker information leakage from the style feature.
Experimental results showed that the proposed model achieves good timbre and
style disentanglement effect, for a given speaker the proposed model could
synthesize speech with any style seen during training, even when the target
style corpus only contains a few hundred training utterances. Furthermore, the
proposed model learns a continuous style representation, which could generate
speech that gradually transits from source style to target style.

## References

  * [1] P. Taylor, _Text-to-speech synthesis_. Cambridge university press, 2009.
  * [2] J. Sotelo, S. Mehri, K. Kumar, J. F. Santos, K. Kastner, A. Courville, and Y. Bengio, “Char2Wav: End-to-end speech synthesis,” in _Proc. ICLR_ , Toulon, 2017.
  * [3] Y. Wang, R. Skerry-Ryan, D. Stanton, Y. Wu, R. J. Weiss, N. Jaitly, Z. Yang, Y. Xiao, Z. Chen, S. Bengio, L. Quoc, and C. R. S. R. Agiomyrgiannakis, Yannis, “Tacotron: Towards end-to-end speech synthesis,” in _Proc. Interspeech_ , Stockholm, 2017.
  * [4] Y. Ren, C. Hu, T. Qin, S. Zhao, Z. Zhao, and T.-Y. Liu, “FastSpeech 2: Fast and High-Quality End-to-End Text to Speech,” _arXiv preprint arXiv:2006.04558_ , 2020.
  * [5] W. Ping, K. Peng, A. Gibiansky, S. O. Arik, A. Kannan, S. Narang, J. Raiman, and J. Miller, “Deep voice 3: Scaling text-to-speech with convolutional sequence learning,” in _arXiv:1710.07654_ , 2017.
  * [6] X. Tan, T. Qin, F. Soong, and T.-Y. Liu, “A Survey on Neural Speech Synthesis,” _arXiv preprint arXiv:2106.15561_ , 2021.
  * [7] X. Tan, J. Chen, H. Liu, J. Cong, C. Zhang, Y. Liu, X. Wang, Y. Leng, Y. Yi, L. He _et al._ , “NaturalSpeech: End-to-End Text to Speech Synthesis with Human-Level Quality,” _arXiv preprint arXiv:2205.04421_ , 2022.
  * [8] S. Pan and L. He, “Cross-speaker style transfer with prosody bottleneck in neural speech synthesis,” _arXiv preprint arXiv:2107.12562_ , 2021.
  * [9] Q. Xie, T. Li, X. Wang, Z. Wang, L. Xie, G. Yu, and G. Wan, “Multi-speaker multi-style text-to-speech synthesis with single-speaker single-style training data scenarios,” _arXiv preprint arXiv:2112.12743_ , 2021.
  * [10] X. An, F. K. Soong, and L. Xie, “Disentangling style and speaker attributes for tts style transfer,” _IEEE/ACM TASLP_ , vol. 30, pp. 646–658, 2022.
  * [11] T. Li, S. Yang, L. Xue, and L. Xie, “Controllable emotion transfer for end-to-end speech synthesis,” in _Proc. ISCSLP 2021_. IEEE, 2021, pp. 1–5.
  * [12] Y. Wang, D. Stanton, Y. Zhang, R.-S. Ryan, E. Battenberg, J. Shor, Y. Xiao, Y. Jia, F. Ren, and R. A. Saurous, “Style tokens: Unsupervised style modeling, control and transfer in end-to-end speech synthesis,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 5180–5189.
  * [13] Y.-J. Zhang, S. Pan, L. He, and Z.-H. Ling, “Learning latent representations for style control and transfer in end-to-end speech synthesis,” in _Proc. ICASSP_. IEEE, 2019, pp. 6945–6949.
  * [14] K. Kumar, R. Kumar, T. de Boissiere, L. Gestin, W. Z. Teoh, J. Sotelo, A. de Brébisson, Y. Bengio, and A. C. Courville, “MelGAN: Generative adversarial networks for conditional waveform synthesis,” in _Proc. NIPS_ , Vancouver, 2019.
  * [15] H. Zen, K. Tokuda, and A. W. Black, “Statistical parametric speech synthesis,” _speech communication_ , vol. 51, no. 11, pp. 1039–1064, 2009\.
  * [16] Q. Xie, X. Tian, G. Liu, K. Song, L. Xie, Z. Wu, H. Li, S. Shi, H. Li, F. Hong, H. Bu, and X. Xu, “The multi-speaker multi-style voice cloning challenge 2021,” in _Proc. ICASSP_. IEEE, 2021\.
  * [17] M. McAuliffe, M. Socolof, S. Mihuc, M. Wagner, and M. Sonderegger, “Montreal Forced Aligner: Trainable Text-Speech Alignment Using Kaldi,” in _Interspeech_ , 2017, pp. 498–502.
  * [18] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _Proc. ICLR_ , San Diego, 2015.
  * [19] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Proc. NIPS_ , Long Beach, 2017.


11institutetext: INAF/Osservatorio Astrofisico di Torino, via Osservatorio 20,
10025 Pino Torinese, Italy 22institutetext: Dipartimento di Fisica, Università
degli Studi di Torino, via Pietro Giuria 1, 10125 Torino, Italy

# The different flavors of extragalactic jets:
Magnetized relativistic flows

P. Rossi 11 G. Bodo 11 S. Massaglia 22 A. Capetti 11

(Received ?? / Accepted ??)

We perform three-dimensional numerical simulations of magnetized relativistic
jets propagating in a uniform density environment in order to study the effect
of the entrainment and the consequent deceleration, extending a previous work
in which magnetic effects were not present. As in previous papers, our aim is
to understand the connection between the jet properties and the resulting
Fanaroff-Riley classification. We consider jets with different low densities,
and therefore low power, and different magnetizations. We find that lower
magnetization jets effectively decelerate to sub-relativistic velocities and
may then result in an FR I morphology on larger scales. At the opposite, in
the higher magnetization cases, the entrainment and consequent deceleration
are substantially reduced.

###### Key Words.:

Magnetohydrodynamics (MHD)– methods: numerical – galaxies: jets

##  1 Introduction

The traditional Fanaroff-Riley classification (Fanaroff & Riley 1974) divides
extragalactic radio-sources in two groups, it is based on morphological
properties, but it is also reflected in radio luminosity. Fanaroff-Riley I (FR
I) are low radio power sources and are characterized by jet-dominated emission
that smoothly extends into the intracluster medium, where large-scale plumes
or tails of diffused emission are observed. The second class, Fanaroff-Riley
II (FR II), is characterized by higher radio power and the sources in this
class show the maximum brightness in the hot spots at the jet termination.
Recently a third class of radio-galaxies has been introduced, representing the
majority of the local radio-loud AGN population. Because they lack the
prominent extended radio structure characteristic of the other Fanaroff-Riley
classes they were called FR 0 (Baldi & Capetti 2009, 2010; Sadler et al. 2014;
Baldi et al. 2015). There is however a minority of FR 0 sources that show
evidences for small-scale jets, usually limited to sizes of a few kpc (Baldi
et al. 2019). This suggests that jets in this class are not even able to
escape from the galaxy core. There are observational evidences that, both in
FR I and in FR II, the jets at their base (at the parsec scale) are
relativistic, with similar Lorentz factors (Giovannini et al. 2001; Celotti &
Ghisellini 2008). The situation is less clear for FR 0 sources, Giovannini et
al. (2023), from VLBI observations, conclude that jets in FR 0 can be, at
most, mildly relativistic, however Massaro et al. (2020) suggest that FR 0
could represent the misaligned counterpart of a significant fraction of BL Lac
objects and therefore they also should be associated to relativistic jets.
This apparent contradiction could be solved by considering a jet with a
transverse velocity structure, with a relativistic spine surrounded by a non
relativistic flow. It is clear that in FR I jets and, possibly also in FR 0
jets, a deceleration to sub-relativistic velocities must occur between the
inner region (at the parsec scale) and the kiloparsec scale, where FR I show
their plume-like morphology (Laing & Bridle 2002a, 2014) and FR 0 disappear.
Massaglia et al. (2016) (hereinafter Mas16) performed high-resolution three-
dimensional simulations of low Mach number jets, showing the onset of
turbulence and the formation of structures very similar to that observed in FR
I sources (see also van der Westhuizen et al. 2019). However, in this study,
the jets are injected already at non-relativistic velocities and one has to
assume that deceleration has already occurred.

Deceleration occurs most likely through the progressive turbulent entrainment
of external material (see e.g. De Young 1993, 2005; Bicknell 1984, 1986, 1994;
Wang et al. 2009), as a consequence of the development of jet instabilities,
in particular Kelvin-Helmholtz (see e.g. Perucho & Martí 2007; Rossi et al.
2008, 2020), but also of other kinds. The role of current-driven instabilities
has been considered, for example, by Tchekhovskoy & Bromberg (2016), Massaglia
et al. (2019, 2022) and Mukherjee et al. (2020). Moreover, it has been
recently pointed out that recollimation shocks may trigger instabilities and
possible entrainment, in particular Gourgouliatos & Komissarov (2018a, b) and
Abolmasov & Bromberg (2023) discussed the development of the relativistic
centrifugal instabilities and Costa et al. (in preparation) consider Kelvin-
Helmholtz instabilities. In addition, we must also mention a second option,
that may be complementary to turbulent entrainment, related to mass loading by
stellar mass loss (Komissarov 1994; Bowman et al. 1996; Laing & Bridle 2002b;
Perucho 2014; Perucho et al. 2014; Perucho 2020; Anglés-Castillo et al. 2021).

In this paper we extend the work by Rossi et al. (2008) and Rossi et al.
(2020) (hereinafter Paper I and Paper II). They performed numerical
simulations of the propagation of relativistic unmagnetized jets, with
different density ratios between the jet and the ambient medium, studying in
detail the entrainment and deceleration process. They found that lighter jets
are more prone to deceleration. In particular, in Paper II, we followed the
jet propagation up to a distance of about 600 initial jet radii, showing how
jets, with a density ratio of $10^{-4}$, are able to reach a quasi-steady
configuration, where the Lorentz factor decreases from 10, at the injection
point, to about 2, close to the jet head. On the contrary, denser jets can
propagate in an almost undisturbed way. The main limitations of the analysis
presented in Papers I and II is the absence of magnetic field, our aim in the
present paper is to overcome this limitation.

In this paper we concentrate our analysis on the role of the density ratio and
of the field strength in determining the jet deceleration process and we will
show that these parameters play a very important role in this process. Of
course, there are also other aspects are likely to play a relevant role on the
jet stability and evolution, such as i) different magnetic field
configurations, ii) different initial Lorentz factors, iii) inhomogeneities in
the ambient medium. However the parameter space is very wide and one has to
proceed in a step-by-step manner.

In Sect. 2 we describe the set of equations to be solved, the numerical
method, the initial and boundary conditions and the parameters of the
simulations. In Sect. 3 we present our results and in Sect. 4 we give
discussion and conclusions.

##  2 Numerical setup

Numerical simulations are carried out by solving the equa- tions of
relativistic ideal magneto-hydrodynamics in 3D. They are expressed by

$\partial_{t}\left(\gamma\rho\right)+\nabla\cdot\left(\gamma\rho\mathbf{v}\right)=0\,,$
(1)

$\partial_{t}\left(\gamma^{2}\rho
h\mathbf{v}+\mathbf{E}\times\mathbf{B}\right)+$
$+\nabla\cdot\left[]\gamma^{2}\rho
h\mathbf{v}\mathbf{v}-\mathbf{E}\mathbf{E}-\mathbf{B}\mathbf{B}+\left(p+u_{em}\right)\mathbf{I}\right]=0\,,$
(2)

$\partial_{t}\left(\gamma^{2}\rho
h-p+u_{em}\right)+\nabla\cdot\left(\gamma^{2}w\mathbf{v}+\mathbf{E}\times\mathbf{B}\right)=0\,,$
(3)

$\partial_{t}\mathbf{B}+\nabla\times\mathbf{E}=0\,,$ (4)

where the electric field $\mathbf{E}$ is provided by the ideal condition
$\mathbf{E}+\mathbf{v}\times\mathbf{B}=0$. In this set of equations, $\rho$ is
the proper density, $\mathbf{v}$ is the velocity three-vector (in units of the
light speed $c$), $\mathbf{B}$ is the laboratory magnetic field (in units of
$\sqrt{4\pi}$), $\gamma=1/\sqrt{1-v^{2}}$ is the Lorentz factor, $h$ is the
relativistic specific enthalpy, $u_{em}=(E^{2}+B^{2})/2$ and $\mathbf{I}$ is
the unit $3\times 3$ tensor. The system of equations (1-4) is completed by
specifying an equation of state relating $h$, $p$ and $\rho$. Following
Mignone et al. (2005), we adopted the Taub-Mathews (TM) equation of state,
with the prescription:

$h=\frac{5}{2}\frac{p}{\rho}+\sqrt{1+\frac{9}{4}\frac{p^{2}}{\rho^{2}}}\,,$
(5)

which closely reproduces the thermodynamics of a Synge gas for a single
species relativistic perfect fluid (notice that we set $c=1$, since the light
speed is our unit of velocity).

In addition to the system (1-4), we also integrate the equation

$\partial_{t}\left(\gamma\rho f\right)+\nabla\cdot\left(\gamma\rho
f\mathbf{v}\right)=0\,,$ (6)

which describes the evolution of a passive tracer $f$, set to unity for the
injected jet material and to zero for the ambient medium. In this way we can
follow the evolution of the spatial distribution of the jet material and its
mixing with the ambient medium.

Simulations are carried out in a Cartesian domain with coordinates in the
range $x\in\left[-L_{x}/2,L_{x}/2\right]$, $y\in\left[0,L_{y}\right]$ and
$z\in\left[-L_{z}/2,L_{z}/2\right]$. At $t=0$ the domain is filled with a
perfect fluid of uniform density and pressure, representing the external
medium initially at rest. In addition, we assume the external medium to be
initially unmagnetized. The assumption of constant density is consistent with
the fact that we focus on the jet propagation inside the galaxy core.

A cylindrical inflow of velocity $v_{j}$, proper density $\rho_{j}$ and a
purely azimuthal magnetic field is constantly fed into the domain at $y=0$,
along the $y$ direction. The jet is characterized by its radius $r_{j}$, the
Lorentz factor $\gamma_{j}$, the ratio $\eta_{j}$ between its proper density
$\rho_{j}$ and the external density and the Mach number defined as
$M_{j}=v_{j}/c_{s0}$, where $c_{s0}$ is the sound speed on the jet axis. The
sound speed $c_{s0}$ is given by the expression valid for the TM equation of
state (Mignone et al. 2005)

$c_{s0}=\frac{p_{0}}{3\rho_{j}h_{0}}\frac{5\rho_{j}h_{0}-p_{0}}{\rho_{j}h_{0}-p_{0}}\,,$
(7)

where $p_{0}$ and $h_{0}$ are respectively the pressure and the specific
enthalpy computed on the jet axis. The specific enthalpy can be derived
through Eq. 5. The injected magnetic field has a configuration corresponding
to a constant current density inside $r=a$ (we take $a=0.75r_{j}$) and zero
outside:

$B_{\varphi}=\cases{-}B_{m}r/a\,,&\text{if}\;r\leq a\\\
-B_{m}a/r\text{if}\;r>a\,,{}$ (8)

where $B_{m}$ is the maximum magnetic field strength reached at $r=a$. The
return current during the jet evolution follows the backflow and the shocked
ambient medium similarly to what is seen in Leismann et al. (2005). Inside the
jet $(r<r_{j}$), we assume total pressure equilibrium, therefore the thermal
pressure inside the jet is given by

$p_{j}(r)=p_{0}-B_{m}^{2}{\rm min}\left(\frac{r^{2}}{a^{2}},1\right)\,.$ (9)

The value of $p_{0}$, as discussed above, is related to the jet Mach number
$M_{j}$, while the value of $B_{m}$, in the jet rest frame, can be obtained
through the magnetization parameter $\sigma_{j}=B_{m}^{2}/\rho_{j}h_{0}$. The
external medium can be characterized the parameter

$\Pi=\frac{p_{a}}{\rho_{a}}\,,$ (10)

where $p_{a}$ and $\rho_{a}$ are respectively the pressure and density of the
external medium.

The equations are non-dimensionalised by using $c$ as the unit of velocity (as
already discussed), the jet radius $r_{j}$ as the unit of length and the
ambient density $\rho_{a}$ as the unit of density. To specify completely the
problem we therefore need the following non-dimensional parameters, introduced
above: the jet Lorentz factor $\gamma_{j}$, the density ratio $\eta_{j}$, the
jet Mach number $M_{j}$, the jet magnetization $\sigma_{j}$, and the non-
dimensional ambient pressure $\Pi$. In our analysis we kept fixed the Lorentz
factor $\gamma_{j}=10$, the Mach number $M_{j}=2.2$ and $\Pi=5\times 10^{-7}$
(corresponding to a temperature of the ambient medium of $5\times 10^{6}$K).
We varied instead the density ratio $\eta_{j}$ and the magnetization
$\sigma_{j}$ for a total of four simulations. All the jets results to be over-
pressured with respect to the ambient medium. The parameters of the
simulations are reported in Table 1, where we also report the pressure
contrast $K$ between the jet and the ambient, the relativistic specific
enthalpy in the jet and the jet power discussed below.

Case | $\eta_{j}$ | $\sigma_{j}$ | $K$ | $h_{0}$ |  $P_{j}$ [erg s-1]
---|---|---|---|---|---
A | $10^{-4}$ | $10^{-2}$ | $20$ | 1.6 | $1.96\times 10^{43}$
B | $10^{-4}$ | $10^{-1}$ | $20$ | 1.6 | $2.04\times 10^{43}$
C | $10^{-5}$ | $10^{-2}$ | $2$ | 1.6 | $1.96\times 10^{42}$
D | $10^{-5}$ | $10^{-1}$ | $2$ | 1.6 | $2.04\times 10^{42}$
Table 1: Parameter set used in the numerical simulation model. The second
column refers to density ratio $\eta_{j}$, the third to the magnetization
$\sigma_{j}$, the fourth gives the derived pressure ratio, the fifth gives the
value of the relativistic specific enthalpy in the jet and the last column to
the jet power evaluated according to Eq. (12)

Numerical integration of the system (1-4) is performed using the ideal
relativistic magnetohydrodynamics module in the PLUTO code (Mignone et al.
2007). For the present application we employed linear reconstruction, second
order Runge-Kutta time stepping and the HLLD Riemann solver (Mignone et al.
2009). Constrained transport (Balsara & Spicer 1999; Londrillo & del Zanna
2004) is used to maintain the condition $\nabla\cdot\mathbf{B}=0$. For the
physical domain extension we choose $L_{x}=300$, $L_{y}=600$ and $L_{z}=300$
and we cover this domain with a grid of $350\times 2400\times 350$ grid
points, which are uniformly spaced in the region around the jet axis
($-15<x<15$, $-15<z<15$) and geometrically stretched outside. The cells are
elongated in the jet directions and we have a resolution of 8 points on the
jet radius in the $x$ and $z$ directions and of 4 points on the jet radius in
the $y$ direction. The boundary conditions are outflow on all the boundaries
except for $y=0$, where, for $r<r_{j}$, we have inflow conditions and, for
$r>r_{j}$ we have reflective conditions.

##  3 Results

Figure 1: Position of the jet head as a function of time, the top panel shows
cases A and B, the bottom panel shows cases C and D.

Figure 2: Three-dimensional isosurfaces of the Lorentz factor at the end of
the simulations. The top-left panel represents case A at $t=23000$, the top-
right panel represents case B at $t=12000$, the bottom-left panel represents
case C at $t=42000$ and, finally, the bottom-right panel represents case D at
$t=48000$. In each panel the top portions displays the entire jet while the
bottom portions displays the last section of the jet between $y=400$ and
$y=600$. The three isosurfaces are for $\gamma=1.1,3,5$.

Our aim is to understand how the deceleration process occurs as a function of
the density contrast and the jet magnetization. We then take two values for
the density contrast, $\eta=10^{-4}$ and $10^{-5}$, and two values for the jet
magnetization, $\sigma_{j}=10^{-2}$ and $10^{-1}$, as shown in Table 1. In
principle the setup of our simulations is scale invariant; however, our
interest is focused on the deceleration that mainly occurs in the first
kiloparsec. This motivates our assumption of constant density of the external
medium because the domain of our simulations can be assumed to be all
contained within the galaxy core. We can then adopt, as in Paper I, a jet
injection radius of the order of one parsec and the external density at
injection of the order of one particle per cubic centimeter (Balmaverde &
Capetti 2006). In the following we will express all quantities in their non-
dimensional values as defined in the previous section. Conversion to physical
units can be obtained with the above assumptions; in particular, the unit of
time is

$\tau\sim 3.26\left(\frac{r_{j}}{1\text{pc}}\right)\text{yrs}\,,$ (11)

and the jet power is

$\begin{split}P_{j}&\sim
h_{0}\left(1+0.5\sigma_{j}\right)\left(\frac{\gamma_{j}}{10}\right)\left(\frac{\eta}{10^{-4}}\right)\\\
&\left(\frac{r_{j}}{1\text{pc}}\right)\left(\frac{n_{0}}{1\text{cm}^{-3}}\right)1.22\times
10^{43}\text{erg s}^{-1}\,,\end{split}$ (12)

where $h_{0}$ is the relativistic specific enthalpy, defined in Eq. 5 and
$n_{0}$ is the external number density. With the assumption of the jet Mach
number $M_{j}=2.2$, the value of $h_{0}$ is 1.6. The jet power for the four
cases is reported in Table 1. We choose the parameters so that the resulting
powers are at about the dividing power between FR I and FR II.

We follow the jet propagation for cases A and B until the bow shock in front
reaches $y=600$, in cases C and D we stop the simulations somewhat earlier
because in these cases the cocoon edge reaches the lateral boundaries. In
these cases the jet head at the end of the simulations reaches $y\sim 550$.

In the following subsection we will discuss the general features of the jet
dynamics, the entrainment and momentun transfer mechanisms and the resulting
details of the jet deceleration.

###  3.1 Jet dynamics

In the two panels of Fig. 1 we show the jet head position as a function of
time for all the cases: more precisely, the top panel shows the cases A and B,
with density ratio $\eta_{j}=10^{-4}$, and the bottom panel shows the cases C
and D, with density ratio $\eta_{j}=10^{-5}$. The dashed black line in each
panel shows the estimated head position obtained by using the theoretical
value of the jet head velocity, in the relativistic case, given by Martí et
al. (1997) (see also Bromberg et al. 2011). This value has been obtained, by
equating the (relativistic) momentum flux of the jet and the ambient medium,
as

$v_{h}=\frac{\gamma_{j}\sqrt{\eta_{R}}}{1+\gamma_{j}\sqrt{\eta_{R}}}\,,\qquad\qquad\eta_{r}=\frac{h_{0}}{h_{a}}\eta\,,$
(13)

where $h_{0}$ and $h_{a}$ are, respectively, the specific relativistic
enthalpies of the jet and of the ambient. In the top panel, we see that case A
(blue curve, lower magnetization $\sigma_{j}=0.01$) follows closely the
theoretical estimate up to $y_{h}\sim 300$ ($t\sim 2500$). Actually the jet-
head velocity reaches values that are slightly higher than the theoretical
estimate and this could be interpreted as the result of jet-head wobbling as
in Aloy et al. (1999) and Perucho et al. (2019). After this initial phase, we
observe a strong deceleration. In Case B (orange curve, higher magnetization
$\sigma_{j}=0.1$), at the beginning, the jet head appears to advance more
slowly than in case A, however the deceleration is much less efficient and the
jet head reaches the end of the domain in about half the time with respect to
case A ($t\sim 12000$ for case B vs $t\sim 23000$ for case B). In the bottom
panel, cases C and D propagate at a much smaller velocity, as it is expected
from the lower density ratio. Case C, with lower magnetization
($\sigma_{j}=0.01$, blue curve), at the beginning, propagates faster than case
D (orange curve, higher magnetization $\sigma_{j}=0.1$), but, at later times,
case C shows a stronger deceleration and the two cases reach the end of the
domain at similar times ($t\sim 45000$).

Figure 3: Longitudinal cut (in the $y-z$ plane) of the distribution of the
Lorentz factor at the final times of the simulations. The four panels refer to
the four different cases (Case A top-left, Case B top-right, Case C bottom-
left, Case D bottom-right). The times are the same as in Fig. 2. The three
coloured bands in the left panels (Cases A and C) highlight the three phases
described in the text. Figure 4: Plot of the entrained mass as a function of
the jet head position for all the cases. The entrained mass is defined as the
mass of the external material moving at $v_{y}>0.1$

Figure 5: Maximum distance from the jet origin at which a given fraction of
momentum flux is carried by material moving at $\gamma\geq 5$ as a function of
time. The four panels refer to the four different cases (Case A top-left; Case
B top-right; Case C bottom-left; Case D bottom-right). In each panel the top
curve shows the jet head position and the other curves represent the fractions
from 10% to 70%. For cases A and C we highlight the areas between subsequent
curves with different colors.

From Fig. 2 we can get an impression of the final structure of the jets. In
the figure, we show, for all the cases, the three-dimensional distribution of
the Lorentz factor, when the jet reached its maximum length. The four
snapshots corresponds to four different times, namely $t\sim 23000$ for case
A, $t\sim 12000$ for case B, $t\sim 42000$ for case C and $t\sim 48000$ for
case D. The differences in times are due to the different jet head propagation
velocities, as we saw above. In the figures, we display three-dimensional
views of the isocontours for three values of the Lorentz gamma factor
($\gamma=1.1,3,5$). Each panel shows both the entire jet and a zoom on the
last 200 length units. The represented domain, in the transverse direction,
extends from $-30$ to $30$, in both directions, for all cases except case A,
that extends from $-40$ to $40$, since, in this case we observe a larger
expansion of the jet. In cases A and C, which are the less magnetized, we have
a behavior similar to the hydrodynamic cases discussed in Paper 2. We observe
a first region, in which the jet seems to propagate straight, then (for
$y\gtrsim 100$ in case A and $y\gtrsim 250$ in case C) we observe the growth
of perturbations, with the appearance of jet wiggling, and eventually the jet
breaks into small high velocity fragments surrounded by slow moving material.
Notice that we do not add any perturbation at the jet inlet and the growth of
perturbations is solely due to the jet-backflow interaction. In case A, the
zoom on the last part of the jet shows that material moving at $\gamma>5$ is
practically absent. The behavior of the more magnetized cases (B and D in the
right panels) is quite different. The high velocity part of the jet remains
more coherent, with some wiggling, but at a larger scale with respect to the
other two cases. The jet fragmentation occurs but the high velocity fragments
have a larger size and are present up to the jet head.

Figure 6: Case A - Top-left panel: maximum pressure found in the computational
domain as a function of time (bottom axis) and the position reached by the jet
head (top axis). The symbols mark the values for each time, while the colored
band identifies the range of variation. Top-right panel: plot of the y
coordinate of maximum pressure, plotted in the top-left panel, as a function
of time (bottom axis) and the position reached by the jet head (top axis).
Bottom-left panel: plot of the maximum Lorentz factor near the jet head as a
function of time (bottom axis) and the position reached by the jet head (top
axis). More precisely the maximum value of $\gamma$ is in the region
$y_{h}-50<y<y_{h}$, where $y_{h}$ is the jet head position. Bottom-right
panel: plot of the jet head velocity as a function of time (bottom axis) and
the position reached by the jet head (top axis).

Figure 7: Same as Fig. 7, but for Case B. The two vertical black lines in the
bottom-right panel indicate one of time interval in which the jet slows down
because of the bending of the jet head as shown in Fig. 10

Figure 8: Same as Fig. 7, but for Case C.

The behavior described above can be better appreciated in Fig. 3, where we
display a longitudinal cut (in the $y-z$ plane) of the distribution of the
Lorentz factor at the final times of the simulations. In Paper II we described
three regions along the jet characterized by a different jet behavior. In the
first region the jet propagates almost undisturbed and we observe the presence
of several recollimation shocks, in the second region perturbations grow
until, in the last region, the jet is broken into a few fast-moving blobs
surrounded by a wide low-velocity envelope. In the present simulations, these
regions are clearly visible for cases A and C: in case A, the unperturbed
region is quite short, up to $y\sim 100$, and characterized by strong
recollimation shocks, the second region spans the jet length from $y\sim 100$
to $y\sim 250$ and for $y\gtrsim 250$ the jet is fragmented; in case C, the
first region ends at $y\sim 250$ and the recollimation shocks are quite weak,
the second region ends at $y\sim 350$ and the fragmentation region is observed
for $y\gtrsim 350$. A comparison with the results of Paper II (see Case C of
Paper II with a density ratio of $10^{-4}$) shows that in case A we have a
much faster growth of perturbations and an earlier jet fragmentation, while
case C behaves similarly to the hydro case. In addition we can also notice
that, in case A, the entrainment region is much wider with respect to case C
(and the hydro cases). Cases B and D, which are more magnetized, show a quite
different behavior, the jet propagates quite straight for larger distances, up
to $y\gtrsim 300$ and then appears to break more suddenly, with the formation
of larger fragments and a much lower entrainment.

###  3.2 Jet entrainment and momentum transfer

The results presented above show that, in cases A and C, the jet progressively
decelerates and the quantity of material moving at high $\gamma$ decreases,
while an increasing amount of material around the jet moves at $\gamma\leq 2$.
Now we can analyze more quantitatively how the entrainment and deceleration
processes occur. In Fig. 4 we show the behavior of the entrained mass as a
function of the jet head position, more precisely we define the entrained mass
as the total mass of the external material moving with $v_{y}>0.1$. The four
curves in the figure refer to the four different cases. Case A, as expected,
is the one with the highest value of the entrained mass. Most of the
entrainment occurs for $400<y_{h}<600$, that corresponds to $7000<t<24000$,
i.e. to the largest fraction of the duration of the jet evolution. For Case C
we had to stop the simulation somewhat earlier with respect to case A because
the lateral bottom part of the cocoon exits the computational grid. Up to this
point, however, case C behaves in a way similar to Case A. For the more
magnetized cases B and D, instead, the entrained mass is much lower, less than
$10\%$ of that of case A.

Figure 9: Same as Fig. 7, but for Case D.  Figure 10: The five panels show
three-dimensional iso-contours of the tracer (for a tracer value of 0.7), for
case B, at five different times in the interval marked by the black vertical
lines in the top-right panel of Fig. 7.

Figure 11: Longitudinal cut (in the $y-z$ plane) of the distribution of the
magnetic field strength at the final times of the simulations. The four panels
refer to the four different cases (Case A top-left, Case B top-right, Case C
bottom-left, Case D bottom-right). The times are the same as in Fig. 2.

In the entrainment process, the jet transfers its momentum to the external
medium. A quantitative view of this process is provided by Fig. 5 where we
plot, as a function of time, the maximum distance at which a given fraction of
the momentum flux is carried by material moving at $\gamma\geq 5$. We start
our discussion by looking at the top left panel, corresponding to case A. In
the figure we see several curves: starting from the top, the blue curve
represents the head position, the region between the top curve and the second
curve refers to a jet section in which the high-velocity material carries a
fraction of momentum flux $<10\%$; the second region refers to a fraction
$10\%-20\%$ and so on; the bottom curve refers to a fraction of $70\%$. Up to
$t\sim 3000$ all the curves follow the behavior of the jet head, subsequently
all the curves tend to flatten and for $t\gtrsim 8000$ they become, on
average, flat. As discussed in Paper II, the flattening shows that the jet has
reached an almost quasi steady state in which it progressively entrains
external material and slows down. We can see that, in this quasi steady state,
the jet up to $y\sim 250$ propagates almost undisturbed, with the high
velocity material carrying most of its momentum, in the region $250\lessapprox
y\lessapprox 400$ the jet progressively transfer its momentum to the entrained
external material, the momentum fraction carried by the high velocity material
drops from $70\%$ to $<10\%$, finally in the region $y>400$ almost all the
initial jet momentum is carried by slow moving material. We can notice that,
as discussed above, the region between $y\sim 250$ and $y\sim 400$ corresponds
to the region where the perturbations progressively increase their amplitude
until the jet becomes fragmented. In the bottom-left panel we show the results
for case C, which has a density contrast of $\eta=10^{-5}$ and the same
magnetization, $\sigma=10^{-2}$ of case A. In this case, a quasi steady state
is almost reached, all the curves tend to flatten, but, on average, they are
still increasing. The process of entrainment and deceleration occurs at a
larger distance with respect to case A, i.e. from $y\sim 350$ to $y\sim 450$
and can be connected, as in case A, with the growth of perturbations and jet
fragmentation that is observed in Fig. 3.

In cases C and D, that have higher magnetization ($\sigma=0.1$), the behavior
is quite different. We already saw that the entrainment process is much less
efficient with respect to the other two cases and therefore the jet does not
lose its momentum almost up to its head. Without entrainment and transfer of
momentum to the external material we would not observe any deceleration and,
in this case, we would expect all the curves to be coincident with the jet
propagation. This is confirmed by the behavior of the different curves in the
right panels of Fig. 5 (top, Case B, $\eta=10^{-4}$; bottom, case D,
$\eta=10^{-5}$), that follow closely the blue line representing the position
of the jet head.

Combining the information presented above, we derive a different scenario for
the low magnetization cases A and C with respect to the high magnetization
case B and D. In cases A and C, the jet is subject to Kelvin-Helmholtz
instability, that drives turbulent entrainment, and progressively releases its
momentum into the ambient medium; it widens its cross-section maintaining a
high gamma value only in a central spine until it breaks and fragments into
high-velocity blobs, surrounded by an envelope at $\gamma\leq 2$. In case A,
however, this process is more efficient and, in the last section of the jet,
the high velocity material is almost absent. In case C, the process is less
efficient it occurs later and fragments of high velocity material are often
present in the last section of the jet, as we will show later. The cases B and
D, with higher magnetization, appear to be less prone to Kelvin-Helmholtz
instabilities, but the jet presents larger scale wiggles that are most likely
due to current-driven instabilities. These differences between low and high
magnetization cases have also been already discussed by Mignone et al. (2010);
Massaglia et al. (2019); Mukherjee et al. (2020) and Massaglia et al. (2022).
In these papers it has been shown how the larger toroidal magnetic fields
stabilize the jet against short wavelength Kelvin-Helmholtz modes that are
responsible of the turbulent entrainment, but lead to current-driven
instabilities that result in larger scale jet wiggling (see also López-
Miralles et al. 2022), in the formation of larger scale fragments, and in a
lower efficiency of entrainment. The development of current-driven instability
is of course enhanced by the presence of only the azimuthal component of
magnetic field. In a more general field configuration, the behavior can be
somewhat different.

###  3.3 Jet deceleration

The entrainment process discussed above may lead to different connected
consequences for the jet dynamics. First, we can have a decrease in the
strength of the terminal shock up to its disappearance. In Massaglia et al.
(2016) the authors interpret this disappearance as a signature of the
transition from FR II to FR I; in fact the presence of a strong terminal shock
can be connected with the presence of hot spots typical of FR II radio
sources. As a second point, the Lorentz factor in the terminal part of the
jet, i.e. in the region close to its head, may tend to decrease, and, finally,
the jet head velocity also decreases. In Figs. 6-9 we give representations of
these different aspects for the four different cases. More precisely, in the
top-left panel of each figure we show a scatter plot of the maximum pressure
in the computational domain, as a function of time (bottom $x$-axis) and jet
head position (top $x$-axis), the coloured band highlights the range of
variation and the general trend. In the top-right panel we plot the $y$
position of the maximum pressure, again as a function of time (bottom
$x$-axis) and jet head position (top $x$-axis). In the bottom-left panel we
plot the maximum value of the Lorentz $\gamma$ near the jet head as a function
of time (bottom $x$-axis) and jet head position (top $x$-axis). More
precisely, in the figure, we plot the maximum value of $\gamma$ found in the
region $y_{h}-50<y<y_{h}$, where $y_{h}$ is the $y$ coordinate of the position
of the jet head. Finally, in the bottom-right panel we show the velocity of
the jet head again as a function of time (lower axis) and of the head position
(upper axis).

In case A (Fig. 6), in the top-left panel, we can observe the presence of two
phases, at first the maximum pressure shows a decreasing trend followed by a
flatter behavior, that starts at $t\gtrsim 10^{4}$, when the jet head has
reached $y\sim 500$. In the flat part, the maximum pressure remains typically
at a low level, with some sporadic spikes. Correspondingly, in the top-right
panel, we observe that, up to $t\sim 10^{4}$ and $y_{h}\sim 500$, the maximum
pressure is found at the jet head, while, for larger times, there are
intervals in which the maximum pressure is observed in the first section of
the jet. This indicates that the terminal shock has disappeared and thus the
jet has acquired a morphology more similar to an FR I. We can notice that the
time at which there is the occurrence of the transition between the two phases
corresponds to the time at which we observe, in the top-left panel, the
transition to a quasi steady-state. The same two phases can also be observed
in the bottom panels. In the bottom-left panel, we can see that, for $t<12000$
($y_{h}<500$), on average, the maximum value of $\gamma$ is around $8$, for
larger times, $t>12000$ ($y_{h}>500$), it drops abruptly to a value of about
2, with some isolated peaks, which are anyway $<4$, except for very few
instances. The jet deceleration is confirmed by the bottom-right panel, where
we observe, in the first phase, strong oscillations of the jet head velocity,
a smoother behavior is instead present in the second phase. As we discussed
earlier, the time of transition between the two phases corresponds to the time
at which the jet reaches an almost quasi steady state. Overall, in case A,
that shows the more pronounced effects of entrainment, the jet head velocity
decreases, from the beginning to the end, by a factor slightly larger than 10.

For case C (Fig. 8), in the top left panel, we similarly observe a two-phase
behavior: in the first phase, up to $t\sim 25000$ and $y_{h}\sim 420$, we
observe a steady decrease, followed again by a flatter part. The range of
variation is however larger than in case A and high pressure peaks occur much
more frequently. Correspondingly, in Fig. 8, top-right panel, we observe that
the intervals in which the terminal shock disappears are very sporadic,
therefore the jet has almost always a FR II-like morphology. Again the
transition between the two phases occurs at the time when the jet has almost
reached a steady state. Similarly, in the bottom-left panel we observe that
for $t<25000$ ($y_{h}<400$) the average value of $\gamma$ is again around $8$,
at larger times it shows strong oscillations with an average value around $4$,
but dropping quite often to values around 2. Also the behavior of the jet head
velocity is similar to case A, with somewhat larger oscillations in the second
phase. In fact, as discussed above, in the quasi steady state, we have more
high velocity fragments that reach the jet head and this leads to the observed
sporadic increases of the jet head velocity. The velocity decrease is somewhat
lower than in Case A.

In the more magnetized cases B and D (Figs. 7 and 9), in the top-left panels,
we do not observe the two-phase behavior, but only a decreasing trend with a
quite large interval of variation and the values of the maximum pressure
remains typically higher than in cases A and C. The position of the maximum
pressure, shown in the top-right panels, is always at the jet head (with only
very few exception for case D) and the maximum values of $\gamma$ stays, for
all times, around $8$, with some sporadic drops at lower values, particularly
for case D. The jet head velocity, in the bottom-right panels, shows a
strongly oscillating behavior for all times. The oscillations can be connected
to the swirlings of the jet head induced by the current-driven kink
instability. This is demonstrated by Fig. 10, where we show a sequence of
three-dimensional views of a tracer iso-contour during the interval of time
marked by the two black vertical lines in the bottom-right panel of Fig. 7
(Case B). We can see that, in correspondence to the decrease of the jet head
velocity, the jet head is deflected at almost $90$ degrees and then the jet
returns to propagate straight and the velocity to increase again. Overall, as
a result of the lower efficiency of the entrainment process, the jet head
velocity for these cases has a lower decrease with respect to the previous
ones.

The structure of the magnetic field appears also quite different between low
and high magnetization cases, as shown in Fig. 11, where we display
longitudinal cuts, in the $y-z$ plane, of the distribution of the magnetic
field intensity. In the two left panels, corresponding to the low
magnetization cases, we can see that for $y\geq 300$ the magnetic field is
characterized by small scale structures. This occurs in the region where
turbulent entrainment is most effective. At the opposite, higher magnetization
cases in the right panels show a much more coherent structure of the magnetic
field concentrated along the jet for almost all its length. This is consistent
with the results of Mukherjee et al. (2020) where similar differences in the
magnetic field distribution are found between the cases dominated by the
Kelvin-Helmholtz instability (as our cases A and C) and the cases dominated by
the current-driven kink instability (as our cases B and D).

##  4 Summary and conclusions

We presented three-dimensional simulations of relativistic magnetized jets
considering four cases with the same Lorentz factor $\gamma=10$, but different
density ratios and different magnetization, for studying the deceleration
process in connection with the FR I-FR II dichotomy. We follow the propagation
up to about 600 pc within the galaxy core, so our assumption of constant
density is well justified. Since in FR I radio sources the jets are
relativistic at their base, but they become sub-relativistic farther out, our
simulations aim, as in Papers I and II, at investigating under which
conditions highly relativistic collimated flows can decelerate to sub-
relativistic velocities within a distance of about one kiloparsec. Here we
extend the analysis performed in Papers I and II, by considering the effects
of magnetic field.

Different kinds of instabilities may play a role in shaping the dynamics
during the jet propagation: in particular Kelvin-Helmholtz instability, as
shown in Papers I and II, can lead to mixing, transfer of momentum and then
jet deceleration. The presence of magnetic field, on the other hand, may have
a strong impact on the instability evolution, in fact it can stabilize the jet
against Kelvin-Helmholtz instabilities, but it can promote the growth of
current-driven instabilities, leading to strong jet bendings and deflections
(see e.g. Mignone et al. 2010).

In Paper I and II we showed how lighter jets, are more prone to deceleration.
In particular the most promising case was Case C of that paper with a density
ratio $\eta=10^{-4}$. In the present paper we then consider cases with
parameters similar to Case C, with $\eta=10^{-4}$, in addition we also
considered cases with a lower density ratio ($\eta=10^{-5}$). For both values
of the density ratios we considered different magnetization. Our results show
striking differences between low and high magnetization cases: while the low
magnetization cases (A and C) show the effects of entrainment and momentum
transfer to the external medium, in high magnetization cases (B and D)
entrainment and momentum transfer appear to be much less efficient. As
discussed above, this difference is due to the magnetic field that has a
stabilizing effect on Kelvin-Helmholtz instabilities and reduces the mixing.
On the other hand, stronger magnetic fields can drive current-driven
instabilities and, in fact, in cases B and D we observe larger scale wigglings
and jet bending. Comparing Cases A and C, we can observe that in case C the
entrainment process is somewhat less efficient than in case A and also it
starts at a more advanced position along the jet. The starting of entrainment
in case A appears to be connected with the presence of strong recollimation
shock at the jet base. The presence of recollimation shocks can favour the
growth of instabilities (Costa et al. in preparation; Abolmasov & Bromberg
2023; Gourgouliatos & Komissarov 2018a, b), in fact, after the shock both the
Lorentz factor and the Mach number of the flow are reduced, making the
conditions more favourable for Kelvin-Helmholtz instabilities. In addition,
one may also observe centrifugal instabilities, recently discussed by
Gourgouliatos & Komissarov (2018a). Case C, at the opposite presents only much
weaker recollimation shocks, this difference may be related to the much higher
jet pressure in Case A. In this context, with respect to the more magnetized
cases B and D, we can note that Mizuno et al. (2015), Martí et al. (2016) and
Moya-Torregrosa et al. (2021) have shown that stronger magnetic field can
avoid the development of strong recollimation shocks, this effect may also
contribute to the increased stability of these cases with respect to KH modes.

In Paper II we showed that low density and therefore low power jets can be
effectively decelerated from highly relativistic to sub-relativistic
velocities within the first kpc. Then they can, on larger scales, display
morphologies typical of FR I radio sources (Massaglia et al. 2016). In this
paper we conclude that the low magnetization cases behave in a way much
similar to the pure hydrodynamic cases discussed in Paper II. More precisely,
these jets reach a quasi steady state in which we observe the progressive
disappearance, along the jet, of material moving at high Lorentz factor. In
the last section of the jet we have a broad flow moving at $\gamma<2$
interspersed with blobs at higher Lorentz factor. In case A, we observe also
the disappearance of the terminal shock, and therefore, from the observational
point of view, of the hot spot. The jet in this case loses the typical
characteristic of a FR II radio source and acquires, already at these early
stages of the jet evolution, an appearance that can be connected to FR I
radio-sources. As discussed by Capetti et al. (2020), the presence of small
scale and/or low brightness jets in FR 0s indicates that FR 0s and FR Is can
be interpreted as two extremes of a continuous population of radio sources,
characterized by a broad distribution of sizes and luminosities of their
extended radio emission. Depending on the physical scale at which the jet
disruption occurs, Case A can be associated with both classes (FR 0 and FR I)
of radio sources with an edge-darkened morphology.

At the opposite, the higher magnetization cases B and D, maintain the jet flow
well collimated up to the end of our simulations. The jets show some
fragmentation, but contrary to the previous cases the fragments are quite
large. In addition, high Lorentz factors are maintained up to the jet head and
we always observe the presence of the terminal shock. Increasing the magnetic
field strength seems to make the transition to an FR I morphology more
difficult even for lower power jets. We have however to consider that we
injected the jets with a pure toroidal field, the presence of a longitudinal
component may somewhat change the behavior. The propagation of jets with a
more complex structure of the magnetic field will be examined in a future
work. Other aspects may as well play a role in the jet entrainment and
deceleration process, one of them is, for example the presence of an
inhomogeneous external medium discussed by Wagner & Bicknell (2011); Wagner et
al. (2012), and Mukherjee et al. (2016), that may lead to a change in the
interaction between the jet and the ambient medium. This will be also examined
in future investigations.

###### Acknowledgements.

This work was supported by PRIN MUR 2022 (grant n. 2022C9TNNX) and by the INAF
Theory Grant Multi scale simulations of relativistic jets. We acknowledge the
INAF-CINECA Accordo Quadro MoU per lo svolgimento di attività congiunta di
ricerca Nuove frontiere in Astrofisica: HPC e Data Exploration di nuova
generazione for the availability of computing resources. We thank M. Perucho
for helpful comments.

## References

  * Abolmasov & Bromberg (2023) Abolmasov, P. & Bromberg, O. 2023, MNRAS, 520, 3009
  * Aloy et al. (1999) Aloy, M. A., Ibáñez, J. M., Martí, J. M., Gómez, J. L., & Müller, E. 1999, ApJ, 523, L125
  * Anglés-Castillo et al. (2021) Anglés-Castillo, A., Perucho, M., Martí, J. M., & Laing, R. A. 2021, MNRAS, 500, 1512
  * Baldi & Capetti (2009) Baldi, R. D. & Capetti, A. 2009, A&A, 508, 603
  * Baldi & Capetti (2010) Baldi, R. D. & Capetti, A. 2010, A&A, 519, A48+
  * Baldi et al. (2015) Baldi, R. D., Capetti, A., & Giovannini, G. 2015, A&A, 576, A38
  * Baldi et al. (2019) Baldi, R. D., Capetti, A., & Giovannini, G. 2019, MNRAS, 482, 2294
  * Balmaverde & Capetti (2006) Balmaverde, B. & Capetti, A. 2006, A&A, 447, 97
  * Balsara & Spicer (1999) Balsara, D. S. & Spicer, D. S. 1999, Journal of Computational Physics, 149, 270
  * Bicknell (1984) Bicknell, G. V. 1984, ApJ, 286, 68
  * Bicknell (1986) Bicknell, G. V. 1986, ApJ, 300, 591
  * Bicknell (1994) Bicknell, G. V. 1994, ApJ, 422, 542
  * Bowman et al. (1996) Bowman, M., Leahy, J. P., & Komissarov, S. S. 1996, MNRAS, 279, 899
  * Bromberg et al. (2011) Bromberg, O., Nakar, E., Piran, T., & Sari, R. 2011, ApJ, 740, 100
  * Capetti et al. (2020) Capetti, A., Brienza, M., Baldi, R. D., et al. 2020, A&A, 642, A107
  * Celotti & Ghisellini (2008) Celotti, A. & Ghisellini, G. 2008, MNRAS, 385, 283
  * Costa et al. ( in preparation) Costa, A., Bodo, G., & Tavecchio, F. in preparation
  * De Young (1993) De Young, D. S. 1993, ApJ, 405, L13
  * De Young (2005) De Young, D. S. 2005, in X-Ray and Radio Connections, ed. L. O. Sjouwerman & K. K. Dyer
  * Fanaroff & Riley (1974) Fanaroff, B. L. & Riley, J. M. 1974, MNRAS, 167, 31P
  * Giovannini et al. (2023) Giovannini, G., Baldi, R. D., Capetti, A., Giroletti, M., & Lico, R. 2023, A&A, 672, A104
  * Giovannini et al. (2001) Giovannini, G., Cotton, W. D., Feretti, L., Lara, L., & Venturi, T. 2001, ApJ, 552, 508
  * Gourgouliatos & Komissarov (2018a) Gourgouliatos, K. N. & Komissarov, S. S. 2018a, Nature Astronomy, 2, 167
  * Gourgouliatos & Komissarov (2018b) Gourgouliatos, K. N. & Komissarov, S. S. 2018b, MNRAS, 475, L125
  * Komissarov (1994) Komissarov, S. S. 1994, MNRAS, 269, 394
  * Laing & Bridle (2002a) Laing, R. A. & Bridle, A. H. 2002a, MNRAS, 336, 1161
  * Laing & Bridle (2002b) Laing, R. A. & Bridle, A. H. 2002b, MNRAS, 336, 328
  * Laing & Bridle (2014) Laing, R. A. & Bridle, A. H. 2014, MNRAS, 437, 3405
  * Leismann et al. (2005) Leismann, T., Antón, L., Aloy, M. A., et al. 2005, A&A, 436, 503
  * Londrillo & del Zanna (2004) Londrillo, P. & del Zanna, L. 2004, Journal of Computational Physics, 195, 17
  * López-Miralles et al. (2022) López-Miralles, J., Perucho, M., Martí, J. M., Migliari, S., & Bosch-Ramon, V. 2022, A&A, 661, A117
  * Martí et al. (1997) Martí, J. M., Müller, E., Font, J. A., Ibáñez, J. M. Z., & Marquina, A. 1997, ApJ, 479, 151
  * Martí et al. (2016) Martí, J. M., Perucho, M., & Gómez, J. L. 2016, ApJ, 831, 163
  * Massaglia et al. (2022) Massaglia, S., Bodo, G., Rossi, P., Capetti, A., & Mignone, A. 2022, A&A, 659, A139
  * Massaglia et al. (2016) Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2016, A&A, 596, A12, Mas16
  * Massaglia et al. (2019) Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2019, A&A, 621, A132
  * Massaro et al. (2020) Massaro, F., Capetti, A., Paggi, A., et al. 2020, ApJ, 900, L34
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
  * Mignone et al. (2005) Mignone, A., Plewa, T., & Bodo, G. 2005, ApJS, 160, 199
  * Mignone et al. (2010) Mignone, A., Rossi, P., Bodo, G., Ferrari, A., & Massaglia, S. 2010, MNRAS, 402, 7
  * Mignone et al. (2009) Mignone, A., Ugliano, M., & Bodo, G. 2009, MNRAS, 393, 1141
  * Mizuno et al. (2015) Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, ApJ, 809, 38
  * Moya-Torregrosa et al. (2021) Moya-Torregrosa, I., Fuentes, A., Martí, J. M., Gómez, J. L., & Perucho, M. 2021, A&A, 650, A60
  * Mukherjee et al. (2016) Mukherjee, D., Bicknell, G. V., Sutherland , R., & Wagner, A. 2016, MNRAS, 461, 967
  * Mukherjee et al. (2020) Mukherjee, D., Bodo, G., Mignone, A., Rossi, P., & Vaidya, B. 2020, MNRAS, 499, 681
  * Perucho (2014) Perucho, M. 2014, in International Journal of Modern Physics Conference Series, Vol. 28, International Journal of Modern Physics Conference Series, 1460165
  * Perucho (2020) Perucho, M. 2020, MNRAS, 494, L22
  * Perucho & Martí (2007) Perucho, M. & Martí, J. M. 2007, MNRAS, 382, 526
  * Perucho et al. (2014) Perucho, M., Martí, J. M., Laing, R. A., & Hardee, P. E. 2014, MNRAS, 441, 1488
  * Perucho et al. (2019) Perucho, M., Martí, J.-M., & Quilis, V. 2019, MNRAS, 482, 3718
  * Rossi et al. (2020) Rossi, P., Bodo, G., Massaglia, S., & Capetti, A. 2020, A&A, 642, A69
  * Rossi et al. (2008) Rossi, P., Mignone, A., Bodo, G., Massaglia, S., & Ferrari, A. 2008, A&A, 488, 795
  * Sadler et al. (2014) Sadler, E. M., Ekers, R. D., Mahony, E. K., Mauch, T., & Murphy, T. 2014, MNRAS, 438, 796
  * Tchekhovskoy & Bromberg (2016) Tchekhovskoy, A. & Bromberg, O. 2016, MNRAS, 461, L46
  * van der Westhuizen et al. (2019) van der Westhuizen, I. P., van Soelen, B., Meintjes, P. J., & Beall, J. H. 2019, MNRAS, 485, 4658
  * Wagner & Bicknell (2011) Wagner, A. Y. & Bicknell, G. V. 2011, ApJ, 728, 29
  * Wagner et al. (2012) Wagner, A. Y., Bicknell, G. V., & Umemura, M. 2012, ApJ, 757, 136
  * Wang et al. (2009) Wang, Y., Kaiser, C. R., Laing, R., et al. 2009, MNRAS, 397, 1113


Automated Counting and Statistical Analysis of Labeled Trees with Degree
Restrictions

Shalosh B. EKHAD and Doron ZEILBERGER

Maple Package

This article is accompanied by a Maple package, ETSIM.txt, available from

https://sites.math.rutgers.edu/~zeilberg/tokhniot/ETSIM.txt .

The web-page of this article,

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/etsim.html ,

contains input and output files, referred to in this paper.

Sneak Preview

In 1887 Arthur Cayley famously proved that there are $n^{n-2}$ labeled trees
on $n$ vertices, but how many trees are there where, for example, all the
vertex-degrees must belong to the set $\\{1,3,4,6,7\\}$? Typing

LTseq($\\{$1,3,4,6,7 $\\}$,50);

in the Maple package ETSIM.txt accompanying the present article, you would get
the first 50 terms immediately. Typing

LTseq($\\{$1,3,4,6,7 $\\}$,1000)[1000];

would give you, just as fast, the $2784$-digit integer that gives you the
exact number of labeled trees with $1000$ vertices such that every vertex
either has one neighbor, or three neighbors, or four neighbors, or six
neighbors, or seven neighbors.

On the other hand, if you want the first $70$ terms of the sequence
enumerating, for example, labeled trees where none of the vertices has degrees
in the set $\\{2,3,5\\}$, just type:

LTseqF($\\{2$,3,5$\\}$,70); .

An Important Formula

Using the general generatingfunctiology for labeled objects, and using
Lagrange Inversion (see [Z] for a lucid and engaging account), one can easily
establish the following theorem.

Important Theorem: Let $P$ be a (finite or infinite) set of positive integers.
Let $a_{P}(n)$ be the number of labeled trees where each vertex must have a
number of neighbors that belongs to $P$, then

$a_{P}(n)$= (n-2)! $\cdot$ Coefficient of $z^{n-2}$ in $\left(\sum_{i\in
P}{z^{i-1}\over(i-1)!}\right)^{n}$ .

Part I: Efficient Counting, Recurrences, and Asymptotics for Enumerating
Labeled trees where all vertices MUST belong to a given Finite Set

In this case the amazing Almkvist-Zeilberger algorithm [AZ] (see [D] for a
lucid and engaging exposition) can easily find a linear recurrence equation
with polynomial coefficients satisfied by $a_{P}(n)$. This is implemented in
procedure LTrec in the Maple package. This also enables one to find
asymptotics. Procedure InfoLT takes care of this. It also gives you the
limiting distribution of the respective participating degrees.

To get an article with information about all possible subsets of
$\\{1,2,\dots,M\\}$ that include $1$ (of course), with the first $K1$ terms of
each sequence displayed, type

LTpaper(M,K1,n); .

Here is the output file when you choose M=7 and K1=30:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oETSIM1a.txt .

The output file

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oETSIM1.txt

only does it for subsets of $\\{1,2,3,4,5\\}$, but also gives the limiting
distribution of the respective participating degrees, that is much more time-
consuming.

Part II: Efficient Counting, and Estimated Asymptotics, for Enumerating
Labeled trees where all vertices must NOT belong to a specified Finite Set

If $F$ is the forbidden set, then in the Important Formula,
$P=\\{1,2,3,\dots\\}\backslash F$, and we have that the number of labeled
trees on $n$ vertices such that none of the vertices have a degree in $F$ is

(n-2)! $\cdot$ Coefficient of $z^{n-2}$ in $\left(e^{z}-\sum_{i\in
F}{z^{i-1}\over(i-1)!}\right)^{n}$ .

If $F$ is a singleton $\\{r\\}$ we have that the the number of labeled trees
on $n$ vertices such that no vertex has $r$ neighbors is

(n-2)! $\cdot$ Coefficient of $z^{n-2}$ in
$\left(e^{z}-{z^{r-1}\over(r-1)!}\right)^{n}$ .

This equals

(n-2)! $\cdot$ Coefficient of $z^{0}$ in
$\left(e^{z}-{z^{r}\over(r-1)!}\right)^{n}\cdot{1\over z^{n-2}}$ ,

that equals

(n-2)! $\cdot$ Coefficient of $z^{0}$ in
$\sum_{k=0}^{n}(-1)^{k}\,{{n}\choose{k}}\,{1\over(r-1)!^{k}}\,{e^{z(n-k)}\over
z^{n-2-r\,k}}$ ,

that equals

$\sum_{k=0}^{\lfloor(n-2)/r\rfloor}\,(-1)^{k}{{n}\choose{k}}\,{1\over(r-1)!^{k}}\,{(n-2)!(n-k)^{n-2-r\,k}\over(n-2-rk)!}\quad,$

that enables a very fast computation. This is implemented in procedure
T1rn(r,n). If the forbidden set has two elements, we have a double sum,
implemented in procedure T2rn(r1,r2,n), and if $F$ has three elements then we
have a triple sum (procedure T3rn(r1,r2,r3,n)).

Procedure LTseqF(F,N) gives the first N terms of the enumerating sequence of
labeled trees avoiding the members of $F$ as vertex-degrees.

Procedure LTFpaper gives information about many cases, together with estimated
(non-rigorous, but very reliable!) asymptotics.

For many such sequences, and estimated asymptotics, see

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oETSIM2.txt .

Part III: Statistical Analysis

Given a labeled tree $T$, Let $X_{d}(T)$ denote the number of vertices that
have degree $d$. It is a random variable, defined on the sample space of all
the $n^{n-2}$ labeled trees on $n$ vertices. In particular $X_{1}(T)$ is the
number of leaves of $T$.

We are interested in the expectation, variance, and higher moments of $X_{d}$.
We will show that for each $d$, $X_{d}$ is asymptotically normal.

We will also explore how $X_{d_{1}}$ and $X_{d_{2}}$ interact, as $n$ goes to
infinity.

These questions are answered by the following very interesting theorem.

Interesting Theorem: Let $d$ be a positive integer, and let $X_{d}(T)$ be the
number of vertices of a labeled tree $T$, that have degree $d$, then we have
the following interesting facts.

$\bullet$

$E[X_{d}]\,=\,{e^{-1}\over(d-1)!}\cdot n+O(1)\quad.$

$\bullet$

$Var(X_{d})=\left({e^{-1}\over(d-1)!}-{(d^{2}-4d+5)e^{-2}\over(d-1)!^{2}}\right)n+O(1)\quad.$

$\bullet$ $X_{d}$ is asymptotically normal.

$\bullet$ If $1<d_{1}<d_{2}$, then

$Cov(X_{d_{1}},X_{d_{2}})\,=\,{2d_{1}+2d_{2}-d_{1}d_{2}-5\over(d_{1}-1)!(d_{2}-1)!e^{2}}\cdot\,n+O(1)\quad.$

$\bullet$ $X_{d_{1}}$ and $X_{d_{2}}$ are joint-asymptotically normal with
(limiting) correlation coefficient

${2d_{1}+2d_{2}-d_{1}d_{2}-5\over\sqrt{(\sqrt{(d_{1}-1)!}-(d_{1}^{2}-4d_{1}+5)e^{-1})\,(\sqrt{(d_{2}-1)!}-(d_{2}^{2}-4d_{2}+5)e^{-1})}}\quad.$

Comment added Feb. 1, 2022: We found out that the expressions for the
expectation and variance of $X_{d}$, as well as the fact that it is (singly-)
asymptotically normal, go back to Alfréd Rény [R], for the $d=1$ case, and to
Amram Meir and John W. Moon [MM] for the general case. See also theorem 7.7
(p. 73) of Moon’s classic monograph [M]. We hope that the result about the
covariance, and the joint asymptotic normality of $X_{d_{1}}$ and $X_{d_{2}}$
are new. We would appreciate any references, of course, in case our hope is
wrong.

Sketch of Proof:

Define the weight of a tree $T$ to be

$wt(T):=\prod_{v\in T}g_{degree(v)}\quad,$

where $g_{1},g_{2},\dots$ are commuting indeterminates. The same reasoning
that lead to the Important Formula tells you that the weight-enumerator of
labeled trees with $n$ vertices is

$(n-2)!\cdot$ Coefficient of $z^{n-2}$ in
$\left(\sum_{i=1}^{\infty}{g_{i}z^{i-1}\over(i-1)!}\right)^{n}$ .

Since we want to focus on $X_{d}$, we set all the $g_{i}$’s to $1$ except for
the active variable, $g_{d}$. This gives us that the probability generating
function (under the uniform distribution on labeled trees) of $X_{d}$, let’s
call it $p_{d}(g_{d})$, is:

${(n-2)!\over n^{n-2}}$ $\cdot$ Coefficient of $z^{n-2}$ in
$\left(e^{z}+(g_{d}-1){z^{d-1}\over(d-1)!}\right)^{n}$ .

The expectation, $E[X_{d}]$, is ${d\over dg_{d}}p(g_{d})\Bigl{|}_{g_{d}=1}$,
and more generally, the $k$-th moment, $E[X_{d}^{k}]$ is $(g_{d}{d\over
dg_{d}})^{k}p(g_{d})\Bigl{|}_{g_{d}=1}$. For any specific $d$ and $k$, Maple
can find explicit expressions, as finite linear combinations of terms of the
form $(n-i)^{n-j}/n^{n-2}$, from which Maple can find, in turn, explicit
expressions for the moments about the mean for each numeric $d$ and each
numeric $k$. From these, for small $k$, one can easily guess explicit
expressions for general $d$, and with a little more effort, one can even get
the computer to do it for symbolic $d$, but still numeric (small) $k$. As $k$
gets larger, the expressions get more and more complicated, but for deducing
limiting behavior, the leading terms suffice, and it is conceivable that one
may be able to do it for both symbolic $d$ and symbolic $k$, but we decided
that we have better things to do.

Similarly for mixed moments. See procedure

LTmom2am(d1,d2,n,k1,k2) ,

that gives you an explicit expression, in $n$, for the mixed $(k_{1},k_{2})$
moment, about the mean, of $X_{d_{1}}$ and $X_{d_{2}}$, but as $k_{1}$ and
$k_{2}$ get larger, the expressions become larger and larger. On the other
hand, Procedure

LTmom2amL(d1,d2,k1,k2) ,

gives us the exact value of the limit, as $n$ goes to infinity, of the former
divided by $n^{\lfloor(d_{1}+d_{2})/2\rfloor}$, from which it is easy to guess
symbolic expressions in $d_{1},d_{2}$ for each specific, numeric, $k_{1}$ and
$k_{2}$.

Since we are experimental mathematicians, doing it for sufficiently many
moments, and comparing the limiting scaled moments to those of the normal
distribution, and the limiting scaled mixed moments to the mixed moments of a
bivariate normal with the above limiting correlation, is good enough for us!

For lots of nice formulas, read the output file

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oETSIM3.txt .

Enjoy!

References

[AZ] Gert Almkvist and Doron Zeilberger, The method of differentiating Under
The integral sign, J. Symbolic Computation 10 (1990), 571-591.
https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/duis.pdf .

[D] Robert Dougherty-Bliss, Integral recurrences from A to Z, Amer. Math.
Monthly, to appear. https://arxiv.org/abs/2102.10170 .

[MM] A. Meir and J.W. Moon, On nodes of degree two in random trees,
Mathematika 15, 188-192.

[M] J.W. Moon, “Counting Labelled Trees”, Candadian Mathematical monographs
#1, Canadian Mathematical Congress, 1970.

[R] A. Rényi, Some remarks on the theory of trees, Publications of the
Mathematical Institute of the Hungarian Academy of Science 4 (1959), 73-85.

[Z] Doron Zeilberger, Lagrange Inversion Without Tears (Analysis) (based on
Henrici), Personal Journal of Shalosh B. Ekhad and Doron Zeilberger, 2002.
https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/lag.html .

Shalosh B. Ekhad and Doron Zeilberger, Department of Mathematics, Rutgers
University (New Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen Rd.,
Piscataway, NJ 08854-8019, USA. Email: [ShaloshBEkhad, DoronZeil] at gmail dot
com .

Exclusively published in the Personal Journal of Shalosh B. Ekhad and Doron
Zeilberger and arxiv.org

First Published: 15 of Shevat, 5782 (Birthday for the Trees; Jan. 17, 2022,
according to their reckoning). This version: Feb. 1, 2022.