variables $z_{i}$ and $z_{i+1}$.

9.16 will be used together with (9.3) to prove 9.3.

######  Proposition 9.16.

Let $u=(u_{1},\ldots,u_{n})$ and $w=(w_{1},\ldots,w_{n})$ be compositions, and
assume that (9.1) or (9.2) holds for $u$ (as appropriate). Then we have the
following.

If $w_{i}>w_{i+1}$ and $u_{i}>u_{i+1}$, we have

(9.4)
$(1-s_{i})(F_{w}^{u}F_{u}(\mathbf{z})+F_{w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}))=\frac{x_{w_{i+1}}(z_{i}-z_{i+1})}{z_{i}(z_{i+1}-x_{w_{i+1}})}(F_{s_{i}w}^{u}F_{u}(\mathbf{z})+F_{s_{i}w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}))$

If $w_{i}>w_{i+1}$ and $u_{i}=u_{i+1}$, we have

(9.5)
$(1-s_{i})F_{w}^{u}F_{u}(\mathbf{z})=\frac{x_{w_{i+1}}(z_{i}-z_{i+1})}{z_{i}(z_{i+1}-x_{w_{i+1}})}(F_{s_{i}w}^{u}F_{u}(\mathbf{z})).$

If $w_{i}=w_{i+1}$ and $u_{i}>u_{i+1}$, we have

(9.6)
$(1-s_{i})(F_{w}^{u}F_{u}(\mathbf{z})+F_{w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}))=0.$

If $w_{i}=w_{i+1}$ and $u_{i}=u_{i+1}$, we have

(9.7) $(1-s_{i})F_{w}^{u}F_{u}(\mathbf{z})=0.$

###### Proof.

We may assume that $max(w)=max(u)+1$ otherwise there is nothing to prove as
$F_{w}^{u}=F_{w}^{s_{i}u}=F_{s_{i}w}^{u}=F_{s_{i}w}^{s_{i}u}=0$. We will prove
(9.4) using 9.12. The proofs of the other three equations are similar using
9.13, 9.14 and 9.15 respectively so we omit them. We divide the proof of (9.4)
into six possible cases.

-Case 1: Suppose $w_{i}<max(w)$ and $u_{i+1}<w_{i+1}$. Then by 9.12 (1), we have $F_{w}^{u}=F_{s_{i}w}^{u}=F_{w}^{s_{i}u}=F_{s_{i}w}^{s_{i}u}=0$ so both sides of (9.4) are zero.

-Case 2: Suppose $w_{i}=max(w)$ and $u_{i+1}<w_{i+1}$. By 9.12 (2), the right-hand side of (9.4) is zero. Using 9.12 (2) and (9.1), the left-hand side of (9.4) is

(9.8) $\displaystyle(1-s_{i})(F_{w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}))$
$\displaystyle=(1-s_{i})\big{(}F^{\prime}(z_{i}-x_{u_{i+1}})z_{i+1}\frac{z_{i}(z_{i+1}-x_{u_{i+1}})}{x_{u_{i+1}}}\frac{(1-s_{i})F_{u}(\mathbf{z})}{z_{i}-z_{i+1}}\big{)}$
$\displaystyle=(1-s_{i})\big{(}F^{\prime}z_{i}z_{i+1}\frac{(z_{i}-x_{u_{i+1}})(z_{i+1}-x_{u_{i+1}})}{x_{u_{i+1}}}\frac{(1-s_{i})F_{u}(\mathbf{z})}{z_{i}-z_{i+1}}\big{)},$

which is zero since we are applying $(1-s_{i})$ to a quantity symmetric in
$z_{i},z_{i+1}$.

-Case 3: Suppose $w_{i}<max(w)$ and $u_{i+1}=w_{i+1}$. By 9.12 (3), (9.4) becomes

$\displaystyle F^{\prime}z_{i}z_{i+1}F_{s_{i}u}(\mathbf{z})$
$\displaystyle=\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}\frac{(1-s_{i})(F^{\prime}z_{i}z_{i+1}F_{u}(\mathbf{z}))}{z_{i}-z_{i+1}},\text{
or equivalently}$ $\displaystyle F_{s_{i}u}(\mathbf{z})$
$\displaystyle=\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}\frac{(1-s_{i})F_{u}(\mathbf{z})}{z_{i}-z_{i+1}},$

which is (9.1) for $u$ and hence true by hypothesis.

-Case 4: Suppose $w_{i}<max(w)$ and $u_{i+1}>w_{i+1}$. The left-hand side of (9.4) is

(9.9)
$(1-s_{i})(F_{w}^{u}F_{u}(\mathbf{z})+F_{w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}))=F^{\prime}z_{i}z_{i+1}(1-s_{i})(F_{u}(\mathbf{z})+F_{s_{i}u}(\mathbf{z})).$

By 9.5, $(1-s_{i})(F_{u}(\mathbf{z})+F_{s_{i}u}(\mathbf{z}))=0$, so the left-
hand side is zero. And by 9.12 (4), the right-hand side of (9.4) is also zero.

-Case 5: Suppose $w_{i}=max(w)$ and $u_{i+1}=w_{i+1}$. Then as shown in (9.8), we have $(1-s_{i})(F_{w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}))=0$. Therefore, multiplying the left-hand side of (9.4) by a constant, and using 9.12 (5), and (9.1), we get

$\displaystyle\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}\frac{(1-s_{i})(F_{w}^{u}F_{u}(\mathbf{z})+F_{w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}))}{z_{i}-z_{i+1}}=\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}\frac{(1-s_{i})(F_{w}^{u}F_{u}(\mathbf{z}))}{z_{i}-z_{i+1}}$
$\displaystyle=\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}\frac{(\frac{(z_{i}-x_{c})z_{i+1}F^{\prime}}{x_{c}}F_{u}(\mathbf{z})-\frac{(z_{i+1}-x_{c})z_{i}F^{\prime}}{x_{c}}s_{i}F_{u}(\mathbf{z}))}{z_{i}-z_{i+1}}$
$\displaystyle=\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}\frac{(\frac{(z_{i+1}-x_{c})z_{i}F^{\prime}}{x_{c}}F_{u}(\mathbf{z})-\frac{(z_{i+1}-x_{c})z_{i}F^{\prime}}{x_{c}}s_{i}F_{u}(\mathbf{z}))}{z_{i}-z_{i+1}}+\frac{(z_{i+1}-x_{w_{i+1}})z_{i}F^{\prime}}{x_{w_{i+1}}}F_{u}(\mathbf{z})$
$\displaystyle=\frac{(z_{i+1}-x_{c})z_{i}F^{\prime}}{x_{c}}\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}\frac{(1-s_{i})F_{u}(\mathbf{z})}{z_{i}-z_{i+1}}+\frac{(z_{i+1}-x_{w_{i+1}})z_{i}F^{\prime}}{x_{w_{i+1}}}F_{u}(\mathbf{z})$
$\displaystyle=\frac{(z_{i+1}-x_{c})z_{i}F^{\prime}}{x_{c}}F_{s_{i}u}(\mathbf{z})+\frac{(z_{i+1}-x_{w_{i+1}})z_{i}F^{\prime}}{x_{w_{i+1}}}F_{u}(\mathbf{z})=F_{s_{i}w}^{u}F_{u}(\mathbf{z})+F_{s_{i}w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}).$

-Case 6: Suppose $w_{i}=max(w)$ and $u_{i+1}>w_{i+1}$. Then by 9.12 (6), we have

$\displaystyle\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}\frac{(1-s_{i})(F_{w}^{u}F_{u}(\mathbf{z})+F_{w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}))}{z_{i}-z_{i+1}}$
$\displaystyle=\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}\frac{(1-s_{i})(\frac{(z_{i}-x_{c})z_{i+1}F^{\prime}}{x_{c}}(F_{u}(\mathbf{z})+F_{s_{i}u}(\mathbf{z})))}{z_{i}-z_{i+1}}$
$\displaystyle=\frac{z_{i}(z_{i+1}-x_{w_{i+1}})}{x_{w_{i+1}}}F^{\prime}(F_{u}(\mathbf{z})+F_{s_{i}u}(\mathbf{z}))=F_{s_{i}w}^{u}F_{u}(\mathbf{z})+F_{s_{i}w}^{s_{i}u}F_{s_{i}u}(\mathbf{z}),$

where we used the fact that $(F_{u}(\mathbf{z})+F_{s_{i}u}(\mathbf{z}))$ is
symmetric in $z_{i}$ and $z_{i+1}$ (9.5). ∎

Now we are ready to prove 9.3.

###### Proof of 9.3.

For a composition $w=(w_{1},\ldots,w_{n})$, it is enough to show that $w$
satisfies (9.1) or (9.2). We use induction on $max(w)$.

Consider the base case $max(w)=2$. There is only one multiline queue
(consisting of one row) of type $w$, so we have

$F_{w}(\mathbf{z})=\prod\limits_{i,w_{i}=1}(z_{i})\prod\limits_{i,w_{i}=2}(z_{i}-x_{1}).$

If $w_{i}=w_{i+1}$, $\eqref{eq:exchange2}$ is immediate. If $w_{i}>w_{i+1}$,
we can write $F_{w}(\mathbf{z})=F^{\prime}\cdot(z_{i}-x_{1})z_{i+1}$ and
$F_{s_{i}w}(\mathbf{z})=F^{\prime}z_{i}(z_{i+1}-x_{1})$ for some $F^{\prime}$
that does not depend on $z_{i}$ and $z_{i+1}$. Therefore the right-hand side
of (9.1) becomes

$\displaystyle\frac{z_{i}(z_{i+1}-x_{1})}{x_{1}}\frac{(1-s_{i})(F^{\prime}(z_{i}-x_{1})z_{i+1})}{z_{i}-z_{i+1}}=\frac{z_{i}(z_{i+1}-x_{1})}{x_{1}}F^{\prime}x_{1}=F^{\prime}z_{i}(z_{i+1}-x_{1})$

which is $F_{s_{i}w}(\mathbf{z})$.

Now consider the case $max(w)>2$ and assume that (9.1) and (9.2) are true for
$u$ such that $max(u)=max(w)-1$. If $w_{i}>w_{i+1}$ then (9.1) follows from
(9.3), (9.4) and (9.5). If $w_{i}=w_{i+1}$ then (9.2) follows from (9.3),
(9.6) and (9.7). ∎

##  10\. The proof of 1.3, the monomial factor conjecture

In this section we will prove 1.3, which gives a formula for the largest
monomial in $\mathbf{x}$ that divides $\psi_{w}$ in the setting where
$y_{i}=0$. We will first give an algebraic argument using the isobaric divided
difference operators to show that
$\eta(w)=\prod_{i=1}^{n-2}x_{i}^{\alpha_{i}(w)+\cdots+\alpha_{n-2}(w)}$ always
divides $\psi_{w}$, see 10.9. We will then use our combinatorial formula for
$\mathbf{z}$-deformed steady state probabilities in terms of multiline queues
(9.3) to show that no greater monomial in $\mathbf{x}$ divides $\psi_{w}$.

######  Definition 10.1.

For a (multivariable) polynomial $p$ and a variable $x$, we write
$x^{d}\mid\mid p$ if $x^{d}\mid p$ and $x^{d+1}\nmid p$, that is, $d$ is the
highest power of $x$ that divides $p$.

Recall from 1.2 that for $w=(w_{1},\ldots,w_{n})\in\operatorname{St}(n)$,
$\alpha_{i}(w)$ is the number of integers greater than $i+1$ among
$\\{i+1=w_{r},w_{r+1},\ldots,w_{s-1},w_{s}=i\\}$, where the subscripts are
taken modulo $n$. To prove 1.3, we will show that

$x_{i}^{\alpha_{i}(w)+\cdots+\alpha_{n-2}(w)}\mid\mid\psi_{w},\hskip
14.22636pt\text{for $1\leq i\leq n-2$}.$

######  Definition 10.2.

Given an integer vector $(b_{1},\ldots,b_{n-2})$ such that $0\leq b_{i}\leq
i$, we define $S(b_{1},\ldots,b_{n-2})\in\operatorname{St}(n)$ as follows. Let
$w^{0}=(1,2,\ldots,n)\in\operatorname{St}(n)$. Then for $1\leq i\leq n-2$, we
recursively construct from $w^{i-1}$ a state $w^{i}$ whose first $n-i-1$
letters are $1,2,\ldots,n-i-1$, by taking the last $b_{i}$ letters of
$w^{i-1}$ and inserting them after $1,2,\ldots,n-i-1$. That is,

$w^{i}=(1,2,\ldots,n-i-1,w^{i-1}_{n-b_{i}+1},\ldots,w^{i-1}_{n-1},w^{i-1}_{n},n-i,w_{n-i+1}^{i-1},\ldots,w_{n-b_{i}}^{i-1}).$

We set $S(b_{1},\ldots,b_{n-2})=w^{n-2}$.

######  Example 10.3.

Using 10.2, we obtain $S(1,1,2,3)=w^{4}\in\operatorname{St}(6)$ as follows.

$\begin{array}[]{lll}w^{0}=(1,2,3,4,5,6)&w^{1}=(1,2,3,4,6,5)&w^{2}=(1,2,3,5,4,6)\\\
w^{3}=(1,2,4,6,3,5)&w^{4}=(1,6,3,5,2,4)&\end{array}$

######  Remark 10.4.

It follows from 10.2 that
$\alpha_{n-1-i}(w^{i})=\alpha_{n-1-i}(w^{i-1})-b_{i}$ and
$\alpha_{j}(w^{i-1})=\alpha_{j}(w^{i})$ if $j\neq n-1-i$.

######  Lemma 10.5.

Let $w=S(b_{1},\ldots,b_{n-2})\in\operatorname{St}(n)$. Then we have

$\alpha_{i}(w)=n-1-i-b_{n-1-i}.$

And for any $w\in\operatorname{St}(n)$, there exists
$(b_{1},\ldots,b_{n-2})\in\mathbb{Z}^{n-2}$ such that $0\leq b_{i}\leq i$ and
$w\sim S(b_{1},\ldots,b_{n-2})$.

###### Proof.

Since $\alpha_{i}(w^{0})=n-1-i$, 10.4 implies that
$\alpha_{i}(w)=n-1-i-b_{n-1-i}.$

There are a total of $(n-1)!$ possible states $S(b_{1},\ldots,b_{n-2})$, which
are all cyclically different since the $\alpha_{i}(w)$’s are cyclically
invariant. So they cover every element in $\operatorname{St}(n)$ up to a
cyclic equivalence. ∎

######  Example 10.6.

Continuing 10.3, if $w:=S(1,1,2,3)=(1,6,3,5,2,4)$ then $\alpha_{1}(w)=1$,
$\alpha_{2}(w)=1$, $\alpha_{3}(w)=1$ and $\alpha_{4}(w)=0$, as claimed in
10.5.

######  Lemma 10.7.

If $\psi_{w^{\prime}}(\mathbf{z})=\pi_{k}(a,b)\psi_{w}(\mathbf{z})$, and if a
monomial $M$ in $x_{1},\ldots,x_{n-1}$ divides $\psi_{w}(\mathbf{z})$, then
$\frac{M}{x_{b}}$ divides $\psi_{w^{\prime}}(\mathbf{z})$.

###### Proof.

This follows from the fact that
$\pi_{k}(a,b)f(\mathbf{z})=\frac{z_{k}(z_{k+1}-x_{b})}{x_{b}}\frac{f(\mathbf{z})-s_{k}f(\mathbf{z})}{z_{k}-z_{k+1}}$,
and each $\psi_{w}(\mathbf{z})$ is a polynomial (4.2). ∎

######  Lemma 10.8.

Consider $w,w^{\prime}\in St(n)$ such that
$w=(1,\ldots,i,w_{i+1},\ldots,w_{n})$ and
$w^{\prime}=(2,\ldots,i,w_{n},w_{i+1},\ldots,w_{n-1},1)$ for $1\leq i\leq
n-1$. Let $w^{\prime\prime}=(1,\ldots i,w_{n},w_{i+1},\ldots,w_{n-1})$. If a
monomial $M$ in $x_{1},\ldots,x_{n-1}$ divides $\psi_{w}(\mathbf{z})$,
$\frac{M}{x_{1}\cdots x_{i}}$ divides both $\psi_{w^{\prime}}(\mathbf{z})$ and
$\psi_{w^{\prime\prime}}(\mathbf{z})$.

###### Proof.

Using 10.7 and the fact that

$\psi_{w^{\prime}}(\mathbf{z})=\pi_{i-1}(w_{n},i)\cdots\pi_{2}(w_{n},3)\pi_{1}(w_{n},2)\pi_{n}(w_{n},1)\psi_{w}(\mathbf{z}),$

we conclude that $\frac{M}{x_{1}\cdots x_{i}}$ divides
$\psi_{w^{\prime}}(\mathbf{z})$. The fact that $\frac{M}{x_{1}\cdots x_{i}}$
divides $\psi_{w^{\prime\prime}}(\mathbf{z})$ follows from 4.2 since
$w^{\prime}$ and $w^{\prime\prime}$ are cyclically equivalent. ∎

######  Proposition 10.9.

For $w\in\operatorname{St}(n)$, we have

$x_{i}^{\alpha_{i}(w)+\cdots+\alpha_{n-2}(w)}\mid\psi_{w},\hskip
14.22636pt\text{for $1\leq i\leq n-2$}.$

###### Proof.

By 4.2, it is enough to show this for any cyclic shift of $w$, so we may
assume $w=S(b_{1},\ldots,b_{n-2})$ for some vector $(b_{1},\ldots,b_{n-2})$.
Using the notation in 10.2, we will use induction on $i$ to show that

(10.1)
$x_{1}^{\alpha_{1}(w^{i})}(x_{1}x_{2})^{\alpha_{2}(w^{i})}\cdots(x_{1}\cdots
x_{n-2})^{\alpha_{n-2}(w^{i})}\mid\psi_{w^{i}}(\mathbf{z})$

Since $w^{0}$ is the identity permutation, we know from the expression for
$\psi_{(1,2,\ldots,n)}(\mathbf{z})$ in 4.2 that
$x_{1}^{n-2}(x_{1}x_{2})^{n-3}\cdots(x_{1}\cdots x_{n-2})^{1}$ divides
$\psi_{w^{0}}(\mathbf{z})$. Now assume (10.1) holds for $i-1$. By repeatedly
applying 10.8, we have that

$\frac{x_{1}^{\alpha_{1}(w^{i-1})}(x_{1}x_{2})^{\alpha_{2}(w^{i-1})}\cdots(x_{1}\cdots
x_{n-2})^{\alpha_{n-2}(w^{i-1})}}{(x_{1}\cdots x_{n-1-i})^{b_{i}}}$

divides $\psi_{{w}^{i}}(\mathbf{z})$. Now 10.4 implies that

$\frac{x_{1}^{\alpha_{1}(w^{i-1})}(x_{1}x_{2})^{\alpha_{2}(w^{i-1})}\cdots(x_{1}\cdots
x_{n-2})^{\alpha_{n-2}(w^{i-1})}}{(x_{1}\cdots
x_{n-1-i})^{b_{i}}}=x_{1}^{\alpha_{1}(w^{i})}(x_{1}x_{2})^{\alpha_{2}(w^{i})}\cdots(x_{1}\cdots
x_{n-2})^{\alpha_{n-2}(w^{i})}.$

Thus we conclude (10.1) is true for $i$. And (10.1) for $i=n-2$ gives the
claim. ∎

Having proved 10.9, our next goal is to show that no larger monomial in
$\mathbf{x}$ divides $\psi_{w}$.

######  Lemma 10.10.

If $x_{i}^{d}\mid\mid\psi_{w}$ then $x_{i}^{d}\mid\mid\psi_{w}(\mathbf{z})$.

###### Proof.

If $x_{i}^{d}\mid\psi_{w}(\mathbf{z})$ then $x_{i}^{d}$ also divides
$\psi_{w}=\operatorname{LC_{z}}(\psi_{w}(\mathbf{z}))$. For the converse,
assume $x_{i}^{d}\mid\psi_{w}$. Then by 8.4,
$x_{i}^{d}\mid\operatorname{wt}(Q)$ for every multiline queue $Q\in MLQ(w)$.
Since $\operatorname{wt}(Q)\mid\operatorname{wt_{z}}(Q)$ we have
$x_{i}^{d}\mid\operatorname{wt_{z}}(Q)$ for every $Q\in MLQ(w)$, which implies
$x_{i}^{d}\mid\psi_{w}(\mathbf{z})$ by 9.3. ∎

######  Lemma 10.11.

Given a multiline queue $Q\in MLQ(w_{1},\ldots,w_{n},1)$, there exists
$\hat{Q}\in MLQ(w_{1}+1,\ldots,w_{n}+1,2,1)$ such that
$x_{1}^{d}\mid\mid\operatorname{wt}(Q)$ if and only if
$x_{1}^{d}\mid\mid\operatorname{wt}(\hat{Q})$.

###### Proof.

Recall the notion of the type of a row of a multiline queue from 8.1. We will
construct $\hat{Q}$ so that if the $r$th row of $Q$ has type
$(u_{1},\ldots,u_{n},1)$ (up to a cyclic shift), the $(r+1)$st row of
$\hat{Q}$ has type $(u_{1}+1,\ldots,u_{n}+1,2,1)$ (up to a cyclic shift).

Let $\hat{Q}^{0}$ be the two-row multiline queue with $n+2$ columns, whose
first and second rows have types $(2,\ldots,2,1,2,2)$ and $(3,\ldots,3,2,1)$.
For $1\leq i\leq n-1$, we inductively construct $\hat{Q}^{i}$ as follows:

  * •

Cyclically shift the columns of $Q$ so that the leftmost entry in the $i$th
row of the resulting multiline queue $Q^{\prime}$ is a ball labeled $1$. Let
$(v_{1},\ldots,v_{k},1,v_{k+1},\ldots,v_{n})$ be the type of the $(i+1)$st row
of $Q^{\prime}$. Add a new $(i+2)$nd row of type
$(v_{1}+1,\ldots,v_{k}+1,2,1,v_{k+1}+1,\ldots,v_{n}+1)$ to $\hat{Q}^{i-1}$,
then cyclically shift the columns of the multiline queue so that the leftmost
entry in the $(i+2)$nd row is a ball labeled $1$. Denote the resulting
multiline queue by $\hat{Q}^{i}$.

We set $\hat{Q}=\hat{Q}^{n-1}$.

$2$$1$$1$$1$$3$$2$$2$$1$$1$$2$$1$$1$$1$$2$$3$$2$$1$$1$$1$$2$$3$$4$$3$$2$$1$$Q$$\hat{Q}^{0}$$\hat{Q}^{1}$$\hat{Q}^{2}=\hat{Q}$
Figure 9. An example of the construction in 10.11 with $n=3$.

By 8.2, the exponent of $x_{1}$ in $\operatorname{wt}(Q)$ is the number of
vacancies in rows $2,\ldots,n$ minus the number of $1$-covered vacancies in
rows $2,\ldots,n$. By construction, the number of vacancies and $1$-covered
vacancies in the $r$th row of $Q$ agrees with the corresponding number of
vacancies in the $(r+1)$st row of $\hat{Q}$ for $r\geq 2$. All the vacancies
in row $2$ of $\hat{Q}$ are $1$-covered, so the exponents of $x_{1}$ in
$\operatorname{wt}(Q)$ and $\operatorname{wt}(\hat{Q})$ are the same. ∎

######  Example 10.12.

Figure 9 shows a multiline queue $Q$ of type $(4,3,2,1)$ as well as
$\hat{Q}^{0}$, $\hat{Q}^{1}$, and $\hat{Q}^{2}=\hat{Q}$, as constructed in the
proof of 10.11. We have

$\displaystyle\operatorname{wt}(Q)=(x_{1})^{1+2}x_{2}^{1}(\frac{x_{2}}{x_{1}})^{2}(\frac{x_{3}}{x_{1}})=x_{2}^{3}x_{3}$
$\displaystyle\operatorname{wt}(\hat{Q})=(x_{1})^{1+2+3}x_{2}^{1+2}x_{3}^{1}(\frac{x_{2}}{x_{1}})^{3}(\frac{x_{3}}{x_{1}})^{2}(\frac{x_{4}}{x_{1}})=x_{2}^{6}x_{3}^{3}x_{4}.$

Note that the exponents of $x_{1}$ in $\operatorname{wt}(Q)$ and
$\operatorname{wt}(\hat{Q})$ are both zero.

######  Proposition 10.13.

For the state $w\in\operatorname{St}(n)$, we have

$x_{1}^{\alpha_{1}(w)+\cdots+\alpha_{n-2}(w)}\mid\mid\psi_{w}.$

###### Proof.

By 10.9, we know that
$x_{1}^{\alpha_{1}(w)+\cdots+\alpha_{n-2}(w)}\mid\psi_{w}.$ So we need to show
that no greater power of $x_{1}$ divides $\psi_{w}$. We use induction, and
suppose this is true for states in $\operatorname{St}(n)$. For the sake of
contradiction, suppose there exists a state $w\in\operatorname{St}(n+1)$ such
that $x_{1}^{\alpha_{1}(w)+\cdots+\alpha_{n-1}(w)+1}\mid\psi_{w}$. By 10.10,
$x_{1}^{\alpha_{1}(w)+\cdots+\alpha_{n-1}(w)+1}\mid\psi_{w}(\mathbf{z})$.
Cyclically shifting the state if necessary, we can assume the last component
of $w$ is a $1$; so write
$w=(w_{1},\ldots,w_{n-k-1},2,w_{n-k},\ldots,w_{n-1},1).$555It is possible that
$w=(w_{1},\ldots,w_{n-1},2,1)$, in which case the following argument
simplifies: there are no $\pi$ operators in the following line, and we take
$w^{\prime}=w$. Then we have

$\psi_{(w_{1},\ldots,w_{n-k-1},2,1,w_{n-k},\ldots,w_{n-1})}(\mathbf{z})=\pi_{n-k+1}(w_{n-k},1)\cdots\pi_{n}(w_{n-1},1)\psi_{w}(\mathbf{z}),$

which implies that
$x_{1}^{\alpha_{1}(w)+\cdots+\alpha_{n-1}(w)+1-k}\mid\psi_{(w_{1},\ldots,w_{n-k-1},2,1,w_{n-k},\ldots,w_{n-1})}(\mathbf{z})$
(by 10.7).

Applying a cyclic shift, we let
$w^{\prime}=(w_{n-k},\ldots,w_{n-1},w_{1},\ldots,w_{n-k-1},2,1);$ we also have
$x_{1}^{\alpha_{1}(w)+\cdots+\alpha_{n-1}(w)+1-k}\mid\psi_{w^{\prime}}$. Since
$\alpha_{1}(w^{\prime})=\alpha_{1}(w)-k$ and
$\alpha_{2}(w^{\prime})=\alpha_{2}(w),\ldots,\alpha_{n-1}(w^{\prime})=\alpha_{n-1}(w)$,
we have

(10.2)
$x_{1}^{a_{1}(w^{\prime})+\cdots+a_{n-1}(w^{\prime})+1}\mid\psi_{w^{\prime}}\text{
for
}w^{\prime}=(w_{n-k},\ldots,w_{n-1},w_{1},\ldots,w_{n-k-1},2,1)\in\operatorname{St}(n+1).$

Consider the state
$w^{\prime\prime}=(w_{n-k}-1,\ldots,w_{n-1}-1,w_{1}-1,\ldots,w_{n-k-1}-1,1)\in\operatorname{St}(n)$.
By the induction hypothesis,
$x_{1}^{\alpha_{1}(w^{\prime\prime})+\cdots+\alpha_{n-2}(w^{\prime\prime})}\mid\mid\psi_{w^{\prime\prime}},$
so there exists a multiline queue $Q$ of type $w^{\prime\prime}$ such that

$x_{1}^{\alpha_{1}(w^{\prime\prime})+\cdots+\alpha_{n-2}(w^{\prime\prime})}\mid\mid\operatorname{wt}(Q).$

By 10.11, there is a multiline queue $Q^{\prime}$ of type $w^{\prime}$ such
that

$x_{1}^{\alpha_{1}(w^{\prime\prime})+\cdots+\alpha_{n-2}(w^{\prime\prime})}\mid\mid\operatorname{wt}(Q^{\prime}).$

Since $\alpha_{1}(w^{\prime})=0$ and
$\alpha_{2}(w^{\prime})=\alpha_{1}(w^{\prime\prime}),\ldots,\alpha_{n-1}(w^{\prime})=\alpha_{n-2}(w^{\prime\prime})$,
we have

$x_{1}^{\alpha_{1}(w^{\prime})+\cdots+\alpha_{n-1}(w^{\prime})}\mid\mid\operatorname{wt}(Q^{\prime}).$

This contradicts (10.2). ∎

######  Lemma 10.14.

For $w\in\operatorname{St}(n)$ and $1\leq i\leq n-3$, let
$w^{\prime}\in\operatorname{St}(n-i)$ be obtained by erasing $1,\ldots,i$ in
$w$ and decreasing the other integers by $i$. Given $Q^{\prime}\in
MLQ(w^{\prime})$, there exists $Q\in MLQ(w)$ such that
$x_{1}^{d}\mid\mid\operatorname{wt}(Q^{\prime})$ if and only if
$x_{i+1}^{d}\mid\mid\operatorname{wt}(Q)$.

###### Proof.

We increase the label of each ball in $Q^{\prime}$ by $i$, and insert trivial
bully paths of type $1,\ldots,i$ in the appropriate positions to get a
multiline queue of type $w$. ∎

######  Example 10.15.

Let $w=(1,4,2,7,6,5,3)\in\operatorname{St}(7)$ and $i=3$. Then
$w^{\prime}=(1,4,3,2)$. Figure 10 shows $Q^{\prime}\in MLQ(w^{\prime})$ on the
left and $Q$ on the right. The weights are

$\displaystyle\operatorname{wt}(Q^{\prime})$
$\displaystyle=x_{1}^{1+2}x_{2}^{1}(\frac{x_{2}}{x_{1}})^{2}(\frac{x_{3}}{x_{1}})=x_{2}^{3}x_{3}$
$\displaystyle\operatorname{wt}(Q)$
$\displaystyle=x_{1}^{1+2+3+4+5}x_{2}^{1+2+3+4}x_{3}^{1+2+3}x_{4}^{1+2}x_{5}^{1}(\frac{x_{5}}{x_{4}})^{2}(\frac{x_{6}}{x_{4}})=x_{1}^{15}x_{2}^{10}x_{3}^{6}x_{5}^{3}x_{6}.$

The exponent of $x_{1}$ in $\operatorname{wt}(Q^{\prime})$ and the exponent of
$x_{4}$ in $\operatorname{wt}(Q)$ are both zero.

$2$$1$$1$$1$$3$$2$$Q^{\prime}$$Q$$1$$2$$4$$6$$5$$3$$1$$2$$3$$5$$4$$1$$2$$3$$4$$1$$3$$2$$1$$2$$1$
Figure 10. An example of 10.14, with $w=(1,4,2,7,6,5,3)$, $i=3$, and
$w^{\prime}=(1,4,3,2)$.

Proof of 1.3. For $w\in\operatorname{St}(n)$ and a number $1\leq i\leq n-3$,
consider the state $w^{\prime}\in\operatorname{St}(n-i)$ obtained by erasing
$1,\ldots,i$ in $w$ and decreasing other integers by $i$. By 10.13,
$x_{1}^{\alpha_{1}(w^{\prime})+\cdots+\alpha_{n-i-2}(w^{\prime})}\mid\mid\psi_{w^{\prime}}$,
hence by 8.4, there exists a multiline queue $Q^{\prime}$ of type $w^{\prime}$
such that

$x_{1}^{\alpha_{1}(w^{\prime})+\cdots+\alpha_{n-i-2}(w^{\prime})}\mid\mid\operatorname{wt}(Q^{\prime}).$

By 10.14, there exists $Q\in MLQ(w)$ such that

(10.3)
$x_{i+1}^{\alpha_{1}(w^{\prime})+\cdots+\alpha_{n-i-2}(w^{\prime})}\mid\mid\operatorname{wt}(Q).$

By the construction of $w^{\prime}$ we have

$\alpha_{i+1}(w)=\alpha_{1}(w^{\prime}),\hskip
11.38092pt\alpha_{i+2}(w)=\alpha_{2}(w^{\prime}),\hskip 11.38092pt\ldots\hskip
11.38092pt\alpha_{n-2}(w)=\alpha_{n-i-2}(w^{\prime}),$

so (10.3) implies

$x_{i+1}^{\alpha_{i+1}(w)+\cdots+\alpha_{n-2}(w)}\mid\mid\operatorname{wt}(Q).$

By 10.9 we have that

$x_{i+1}^{\alpha_{i+1}(w)+\cdots+\alpha_{n-2}(w)}\mid\mid\psi_{w}$

for $1\leq i\leq n-3$. Combining this with 10.13, we conclude that the largest
monomial that can be factored out of $\psi_{w}$ is
$\prod\limits_{i=1}^{n-2}x_{i}^{\alpha_{i}(w)+\cdots+\alpha_{n-2}(w)}$, as
desired. Recall from 10.5 that if two states
$w,w^{\prime}\in\operatorname{St}(n)$ have the same
$\alpha_{i}(w)=\alpha_{i}(w^{\prime})$ for all $1\leq i\leq n-2$ then $w\sim
w^{\prime}$. This completes the proof. ∎

##  11\. Future questions

There are various natural questions that arise from this work.

######  Problem 11.1.

Find a combinatorial formula for steady state probabilities of the
inhomogeneous TASEP (e.g. using multiline queues) in the case where $y_{i}\neq
0$.

######  Problem 11.2.

Find a purely combinatorial proof of 1.3.

######  Problem 11.3.

Prove that for any permutation $w$, the steady state probability $\psi_{w}$ in
the inhomogeneous TASEP can be written as a positive sum of double Schubert
polynomials $\mathfrak{S}_{w}(x_{1},\ldots,x_{n};y_{1},\ldots,y_{n})$
multiplied by some linear factors $(x_{i}-y_{j})$.

11.3 is open even when $y_{i}=0$ for all $i$. One way to approach this problem
would be to find a geometric interpretation of each steady state probability.
Recall that Schubert polynomials represent cohomology classes of Schubert
varieties in the complete flag variety. Can one associate a variety to each
state whose cohomology class is represented by the corresponding
(unnormalized) steady state probability?

##  Appendix A Technical results for the proof of 5.9

In this section, we collect some technical results (A.4 and A.5) which are
used for the proof of 5.9.

######  Lemma A.1.

[Can16, Lemma 22] Let $K(x_{1};x_{2},\ldots,x_{m+1})$ be a rational expression
in $x_{1},\ldots,x_{m+1}$ which is symmetric in the variables
$x_{2},\ldots,x_{m+1}$. We have

$\partial_{m}\cdots\partial_{1}K=\sum_{i=1}^{m+1}\frac{K(x_{i};x_{1},\ldots,\hat{x_{i}},\ldots,x_{m+1})}{\prod\limits_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{m+1}(x_{i}-x_{j})}.$

Now we rewrite (5.1). In what follows, we let $[m]$ denote
$\\{1,2,\ldots,m\\}$.

######  Lemma A.2.

For $\lambda\in\operatorname{Val}(n)$, let $\operatorname{mul}(\lambda)\geq k$
for some $k$. Let $\tilde{\lambda}$ denote the partition obtained by deleting
the first $k$ parts of $\lambda$. Then
$\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})$ is equal to

$\sum\limits_{\begin{subarray}{c}I\subseteq[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\\\
|I|=k\end{subarray}}\frac{\displaystyle{\mathfrak{S}_{\tilde{\lambda}}^{n-k}(\sigma^{\lambda_{1}-\lambda_{k+1}+k}(\mathbf{z});\mathbf{x}_{\hat{I}};\mathbf{y})\prod\limits_{\begin{subarray}{c}1\leq
l\leq n-\operatorname{mul}(\lambda)\\\ i\in
I\end{subarray}}(x_{i}-y_{l})\prod\limits_{\begin{subarray}{c}1\leq
i\leq\lambda_{1}-\lambda_{k+1}+k\\\
j\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I\end{subarray}}(z_{i}-x_{j})}}{\displaystyle{\prod\limits_{\begin{subarray}{c}i\in
I\\\ j\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I\end{subarray}}(x_{i}-x_{j})}}$

###### Proof.

Note that the expression in A.2 implies that
$\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})$ is symmetric in
variables $x_{1},\ldots,x_{n-\lambda_{1}-\operatorname{mul}(\lambda)+k}$.
Since we can take $k=\operatorname{mul}(\lambda)$, A.2 implies that
$\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})$ is symmetric in
variables $x_{1},\ldots,x_{n-\lambda_{1}}$. Now we use induction on $n$.
Suppose the statement holds for partitions in $\operatorname{Val}(n-1)$.

Let $\lambda^{\prime}$ to be the partition obtained by deleting the first part
of $\lambda$. By the induction hypothesis,
$\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\mathbf{z};\mathbf{x};\mathbf{y})$ is
symmetric in variables $x_{1},\ldots,x_{(n-1)-\lambda^{\prime}_{1}}.$ Note
that we have $(n-1-\lambda^{\prime}_{1})\geq
n-\lambda_{1}-\operatorname{mul}(\lambda)$. Thus we rewrite (5.1) using A.1 as
follows

(A.1)
$\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})=\sum\limits_{i=1}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}\frac{\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\sigma^{\lambda_{1}-\lambda_{2}+1}(\mathbf{z});\mathbf{x}_{\hat{i}};\mathbf{y})\prod\limits_{l=1}^{n-\operatorname{mul}(\lambda)}(x_{i}-y_{l})\prod\limits_{l=1}^{(\lambda_{1}-\lambda_{2}+1)}\prod\limits_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}(z_{l}-x_{j})}{\prod\limits_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}(x_{i}-x_{j})},$

which implies the $k=1$ case of A.2.

Assume $k>1$. In this case we have $\operatorname{mul}(\lambda)>1$, which
implies $\operatorname{mul}(\lambda^{\prime})=\operatorname{mul}(\lambda)-1$
and $\lambda^{\prime}_{1}=\lambda_{2}=\lambda_{1}$. Applying A.2 to
$\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\mathbf{z};\mathbf{x};\mathbf{y})$ and
$(k-1)$ gives

(A.2)
$\displaystyle\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\mathbf{z};\mathbf{x};\mathbf{y})=$
$\displaystyle\sum\limits_{\begin{subarray}{c}I\subseteq[n-1-\lambda_{1}-\operatorname{mul}(\lambda)+k]\\\
|I|=k-1\end{subarray}}\frac{\mathfrak{S}_{\tilde{\lambda}}^{n-k}(\sigma^{\lambda_{1}-\lambda_{k+1}+k-1}(\mathbf{z});\mathbf{x}_{\hat{I}};\mathbf{y})\prod\limits_{\begin{subarray}{c}1\leq
l\leq n-\operatorname{mul}(\lambda)\\\ i\in
I\end{subarray}}(x_{i}-y_{l})\prod\limits_{\begin{subarray}{c}1\leq
i\leq\lambda_{1}-\lambda_{k+1}+k-1\\\
j\in[n-1-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I\end{subarray}}(z_{i}-x_{j})}{\prod\limits_{\begin{subarray}{c}i\in I\\\
j\in[n-1-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I\end{subarray}}(x_{i}-x_{j})},$

and rewriting (A.1) using the fact $\lambda_{1}=\lambda_{2}$ gives

(A.3)
$\displaystyle\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})=\sum\limits_{j=1}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}\frac{\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\sigma(\mathbf{z});\mathbf{x}_{\hat{j}};\mathbf{y})\prod\limits_{l=1}^{n-\operatorname{mul}(\lambda)}(x_{j}-y_{l})\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq
j\end{subarray}}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}(z_{1}-x_{m})}{\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq
j\end{subarray}}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}(x_{j}-x_{m})}.$

Plugging (A.2) into (A.3) gives

(A.4)
$\displaystyle\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})=\sum\limits_{\begin{subarray}{c}I\subseteq[n-1-\lambda_{1}-\operatorname{mul}(\lambda)+k]\\\
|I|=k-1\end{subarray}}\sum\limits_{j=1}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}f_{I}(\sigma^{\lambda_{1}-\lambda_{k+1}+k}(\mathbf{z});\mathbf{x}_{\hat{j}};\mathbf{y})M_{I}(\mathbf{z};\mathbf{x}_{\hat{j}};\mathbf{y})\prod\limits_{l=1}^{n-\operatorname{mul}(\lambda)}(x_{j}-y_{l})$

where

$\displaystyle f_{I}(\mathbf{z};\mathbf{x};\mathbf{y})$
$\displaystyle=\mathfrak{S}_{\tilde{\lambda}}^{n-k}(\mathbf{z};\mathbf{x}_{\hat{I}};\mathbf{y})$
$\displaystyle M_{I}(\mathbf{z};\mathbf{x};\mathbf{y})$
$\displaystyle=\frac{\prod\limits_{\begin{subarray}{c}1\leq l\leq
n-\operatorname{mul}(\lambda)\\\ i\in
I\end{subarray}}(x_{i}-y_{l})\prod\limits_{\begin{subarray}{c}2\leq
i\leq\lambda_{1}-\lambda_{k+1}+k\\\
m\in[n-1-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I\end{subarray}}(z_{i}-x_{m})\prod\limits_{m=1}^{n-\lambda_{1}-\operatorname{mul}(\lambda)}(z_{1}-x_{m})}{\prod\limits_{\begin{subarray}{c}i\in
I\\\ m\in[n-1-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I\end{subarray}}(x_{i}-x_{m})\prod\limits_{m=1}^{n-\lambda_{1}-\operatorname{mul}(\lambda)}(x_{j}-x_{m})}.$

For a fixed
$I_{0}=\\{i_{1}<\cdots<i_{k}\\}\subseteq[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]$,
to have

$f_{I}(\sigma^{\lambda_{1}-\lambda_{k+1}+k}(\mathbf{z});\mathbf{x}_{\hat{j}};\mathbf{y})=\mathfrak{S}_{\tilde{\lambda}}^{n-k}(\sigma^{\lambda_{1}-\lambda_{k+1}+k}(\mathbf{z});\mathbf{x}_{\hat{I_{0}}};\mathbf{y})$

we need to take $I=I_{h}=\\{i_{1}<\cdots<i_{h-1}<i_{h+1}-1<\cdots<i_{k}-1\\}$
and $j=i_{h}\leq n-\lambda_{1}-\operatorname{mul}(\lambda)+1$. So taking the
coefficient of
$\mathfrak{S}_{\tilde{\lambda}}^{n-k}(\sigma^{\lambda_{1}-\lambda_{k+1}+k}(\mathbf{z});\mathbf{x}_{\hat{I_{0}}};\mathbf{y})$
in (A.4) gives

$\sum_{\begin{subarray}{c}h\geq 1\\\ i_{h}\leq
n-\lambda_{1}-\operatorname{mul}(\lambda)+1\end{subarray}}M_{I_{h}}(\mathbf{z};\mathbf{x}_{\hat{i_{h}}};\mathbf{y})\prod\limits_{l=1}^{n-\operatorname{mul}(\lambda)}(x_{i_{h}}-y_{l}).$

Note that we have

$\displaystyle
M_{I_{h}}(\mathbf{z};\mathbf{x}_{\hat{i_{h}}};\mathbf{y})\prod\limits_{l=1}^{n-\operatorname{mul}(\lambda)}(x_{i_{h}}-y_{l})$
$\displaystyle=\frac{\prod\limits_{\begin{subarray}{c}1\leq l\leq
n-\operatorname{mul}(\lambda)\\\ i\in
I_{0}\end{subarray}}(x_{i}-y_{l})\prod\limits_{\begin{subarray}{c}2\leq
i\leq\lambda_{1}-\lambda_{k+1}+k\\\
m\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I_{0}\end{subarray}}(z_{i}-x_{m})\prod\limits_{\begin{subarray}{c}m=1\\\ m\neq
i_{h}\end{subarray}}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}(z_{1}-x_{m})}{\prod\limits_{\begin{subarray}{c}i\in
I_{0},i\neq i_{h}\\\
m\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I_{0}\end{subarray}}(x_{i}-x_{m})\prod\limits_{\begin{subarray}{c}m=1\\\ m\neq
i_{h}\end{subarray}}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}(x_{i_{h}}-x_{m})}.$

We claim

(A.5) $\frac{\prod\limits_{\begin{subarray}{c}1\leq l\leq
n-\operatorname{mul}(\lambda)\\\ i\in
I_{0}\end{subarray}}(x_{i}-y_{l})\prod\limits_{\begin{subarray}{c}1\leq
i\leq\lambda_{1}-\lambda_{k+1}+k\\\
m\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I_{0}\end{subarray}}(z_{i}-x_{m})}{\prod\limits_{\begin{subarray}{c}i\in
I_{0}\\\ m\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I_{0}\end{subarray}}(x_{i}-x_{m})}=\sum_{\begin{subarray}{c}h\geq 1\\\
i_{h}\leq
n-\lambda_{1}-\operatorname{mul}(\lambda)+1\end{subarray}}M_{I_{h}}(\mathbf{z};\mathbf{x}_{\hat{i_{h}}};\mathbf{y})\prod\limits_{l=1}^{n-\operatorname{mul}(\lambda)}(x_{i_{h}}-y_{l}).$

We view both sides as polynomials in $z_{1}$ of degree at most
$(n-\lambda_{1}-\operatorname{mul}(\lambda))$. So it is enough to show that
they coincide when we plug in
$z_{1}=x_{1},\ldots,x_{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}$. If we
plug in $z=x_{m}$ for $m\notin I_{0}$ then both sides are zero. If we plug in
$z_{1}=x_{i_{t}}$ for some $i_{t}\in I_{0}$ then only $h=t$ on the right hand
side does not vanish and we have

$\big{(}M_{I_{t}}(\mathbf{z};\mathbf{x}_{\hat{i_{t}}};\mathbf{y})\prod\limits_{l=1}^{n-\operatorname{mul}(\lambda)}(x_{i_{t}}-y_{l})\big{)}|_{z_{1}=x_{i_{t}}}=\frac{\prod\limits_{\begin{subarray}{c}1\leq
l\leq n-\operatorname{mul}(\lambda)\\\ i\in
I_{0}\end{subarray}}(x_{i}-y_{l})\prod\limits_{\begin{subarray}{c}2\leq
i\leq\lambda_{1}-\lambda_{k+1}+k\\\
m\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I_{0}\end{subarray}}(z_{i}-x_{m})}{\prod\limits_{\begin{subarray}{c}i\in
I_{0},i\neq i_{t}\\\
m\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I_{0}\end{subarray}}(x_{i}-x_{m})}.$

This is the same as the left hand side of (A.5) evaluated at
$z_{1}=x_{i_{t}}$.

In conclusion (A.4) becomes

$\displaystyle\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})=$
$\displaystyle\sum\limits_{\begin{subarray}{c}I_{0}\subseteq[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\\\
|I_{0}|=k\end{subarray}}\frac{\mathfrak{S}_{\tilde{\lambda}}^{n-k}(\sigma^{\lambda_{1}-\lambda_{k+1}+k}(\mathbf{z});\mathbf{x}_{\hat{I_{0}}};\mathbf{y})\prod\limits_{\begin{subarray}{c}1\leq
l\leq n-\operatorname{mul}(\lambda)\\\ i\in
I_{0}\end{subarray}}(x_{i}-y_{l})\prod\limits_{\begin{subarray}{c}1\leq
i\leq\lambda_{1}-\lambda_{k+1}+k\\\
m\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I_{0}\end{subarray}}(z_{i}-x_{m})}{\prod\limits_{\begin{subarray}{c}i\in
I_{0}\\\ m\in[n-\lambda_{1}-\operatorname{mul}(\lambda)+k]\setminus
I_{0}\end{subarray}}(x_{i}-x_{m})},$

which completes the proof. ∎

######  Proposition A.3.

Choose $\lambda\in\operatorname{Val}(n)$ and $a$ such that $1\leq a\leq
n-\lambda_{1}$. If we specialize $z_{1}:=x_{a}$ in the $\mathbf{z}$-Schubert
polynomial $\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})$, we
obtain

(A.6)
$\displaystyle\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y}))|_{z_{1}=x_{a}}=$
$\displaystyle\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\sigma^{\lambda_{1}-\lambda_{2}+1}(\mathbf{z});\mathbf{x}_{\hat{a}};\mathbf{y}))\prod\limits_{l=1}^{n-\operatorname{mul}(\lambda)}(x_{a}-y_{l})\prod\limits_{i=2}^{(\lambda_{1}-\lambda_{2}+1)}\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq
a\end{subarray}}^{n-\lambda_{1}-\operatorname{mul}(\lambda)+1}(z_{i}-x_{m}),$

where $\lambda^{\prime}$ is the partition obtained by deleting the first part
of $\lambda$. If $\operatorname{length}(\lambda)=1$ then we regard
$\lambda_{2}=0$.

###### Proof.

Taking $k=\operatorname{mul}(\lambda)$ in (A.2), we see that
$\mathfrak{S}_{\lambda}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})$ is symmetric in
variables $x_{1},\ldots,x_{n-\lambda_{1}}$. So it is enough to show (A.6) for
$a=1$, which follows from (A.1). ∎

Note that in A.4 and A.5, the operator $\partial_{i}^{\mathbf{z}}$ _acts on
$\mathbf{z}$-variables_.

######  Proposition A.4.

Let $\lambda\in\operatorname{Val}(n)$ such that $\lambda_{1}>\lambda_{2}$, and
write $\lambda=(\lambda_{1},\lambda^{\prime})$ for some $\lambda^{\prime}$.
Then

(A.7)
$\displaystyle\mathfrak{S}_{(\lambda_{1},\lambda^{\prime})}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})\prod\limits_{i=3}^{\lambda_{1}-\lambda_{2}+1}(z_{i}-x_{n+1-\lambda_{1}})$
$\displaystyle=\partial^{\mathbf{z}}_{1}\big{(}\mathfrak{S}_{(\lambda_{1}-1,\lambda^{\prime})}^{n}(\sigma(\mathbf{z});\mathbf{x};\mathbf{y})\prod\limits_{i=1}^{n-\lambda_{1}}(z_{1}-x_{i})\prod\limits_{i=1}^{\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))-1}(z_{2}-y_{n-i})\big{)}.$

###### Proof.

We will refer to the left and right-hand sides of (A.7) as LHS and RHS. To
prove they are equal, we will view each side as a polynomial in $z_{1}$ and
analyze its degree to use interpolation. Note that
$\mathfrak{S}_{(\lambda_{1}-1,\lambda^{\prime})}^{n}(\sigma(\mathbf{z});\mathbf{x};\mathbf{y})$
does not depend on $z_{1}$ so exchanging variables $z_{1}$ and $z_{2}$ is the
same as plugging in $z_{1}$ in place of $z_{2}$. So we get

$\displaystyle RHS=\frac{1}{z_{1}-z_{2}}$
$\displaystyle\big{[}\mathfrak{S}_{(\lambda_{1}-1,\lambda^{\prime})}^{n}(\sigma(\mathbf{z});\mathbf{x};\mathbf{y})\prod\limits_{i=1}^{n-\lambda_{1}}(z_{1}-x_{i})\prod\limits_{i=1}^{\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))-1}(z_{2}-y_{n-i})$
$\displaystyle-\mathfrak{S}_{(\lambda_{1}-1,\lambda^{\prime})}^{n}(\sigma(\mathbf{z});\mathbf{x};\mathbf{y})|_{z_{2}=z_{1}}\prod\limits_{i=1}^{n-\lambda_{1}}(z_{2}-x_{i})\prod\limits_{i=1}^{\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))-1}(z_{1}-y_{n-i})\big{]}.$

By (5.1),
$\mathfrak{S}_{(\lambda_{1}-1,\lambda)}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})$
has degree at most
$(n-(\lambda_{1}-1)-\operatorname{mul}((\lambda-1,\lambda^{\prime})))$ in
$z_{1}$. So the numerator of $RHS$ has degree at most $(n-\lambda_{1})$ in
$z_{1}$, which implies that $RHS$ has degree at most $(n-\lambda_{1}-1)$ in
$z_{1}$. Meanwhile, LHS has degree at most
$(n-\lambda_{1}-\operatorname{mul}(\lambda))=(n-\lambda_{1}-1)$ in $z_{1}$ by
(5.1). So it is enough to show that LHS and RHS coincide when we plug in
$z_{1}=x_{h}$ for any $1\leq h\leq n-\lambda_{1}$. We have

$\displaystyle RHS|_{z_{1}=x_{h}}$
$\displaystyle=\frac{-\mathfrak{S}_{(\lambda_{1}-1,\lambda^{\prime})}^{n}(\sigma(\mathbf{z});\mathbf{x};\mathbf{y})|_{z_{2}=x_{h}}\prod\limits_{i=1}^{n-\lambda_{1}}(z_{2}-x_{i})\prod\limits_{i=1}^{\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))-1}(x_{h}-y_{n-i})}{x_{h}-z_{2}}$
$\displaystyle=\mathfrak{S}_{(\lambda_{1}-1,\lambda^{\prime})}^{n}(\sigma(\mathbf{z});\mathbf{x};\mathbf{y})|_{z_{2}=x_{h}}\prod\limits_{\begin{subarray}{c}i=1\\\
i\neq
h\end{subarray}}^{n-\lambda_{1}}(z_{2}-x_{i})\prod\limits_{i=1}^{\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))-1}(x_{h}-y_{n-i}).$

And by A.3 we have

$\displaystyle\mathfrak{S}_{(\lambda_{1}-1,\lambda^{\prime})}^{n}(\mathbf{z};\mathbf{x};\mathbf{y}))|_{z_{1}=x_{h}}$
$\displaystyle=\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\sigma^{\lambda_{1}-\lambda_{2}}(\mathbf{z});\mathbf{x}_{\hat{h}};\mathbf{y}))\prod\limits_{l=1}^{n-\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))}(x_{h}-y_{l})\prod\limits_{i=2}^{(\lambda_{1}-\lambda_{2})}\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq
h\end{subarray}}^{n-\lambda_{1}-\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))+2}(z_{i}-x_{m})$
$\displaystyle=\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\sigma^{\lambda_{1}-\lambda_{2}}(\mathbf{z});\mathbf{x}_{\hat{h}};\mathbf{y}))\prod\limits_{l=1}^{n-\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))}(x_{h}-y_{l})\prod\limits_{i=2}^{(\lambda_{1}-\lambda_{2})}\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq h\end{subarray}}^{n-\lambda_{1}+1}(z_{i}-x_{m}),$

where the last equality uses the fact that when
$\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))>1$ we have
$\lambda_{1}-\lambda_{2}=1$ so the product over $i=2$ to
$\lambda_{1}-\lambda_{2}$ is vacuous. Shifting variables by
$\mathbf{z}\rightarrow\sigma(\mathbf{z})$ we deduce

$\displaystyle\mathfrak{S}_{(\lambda_{1}-1,\lambda)}^{n}(\sigma(\mathbf{z});\mathbf{x};\mathbf{y})|_{z_{2}=x_{h}}=\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\sigma^{\lambda_{1}-\lambda_{2}+1}(\mathbf{z});\mathbf{x}_{\hat{h}};\mathbf{y})\prod\limits_{l=1}^{n-\operatorname{mul}((\lambda_{1}-1,\lambda^{\prime}))}(x_{h}-y_{l})\prod\limits_{i=3}^{(\lambda_{1}-\lambda_{2}+1)}\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq h\end{subarray}}^{n-\lambda_{1}-1}(z_{i}-x_{m}).$

Plugging this into $RHS|_{z_{1}=x_{h}}$ gives

$\displaystyle RHS|_{z_{1}=x_{h}}$
$\displaystyle=\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\sigma^{\lambda_{1}-\lambda_{2}+1}(\mathbf{z});\mathbf{x}_{\hat{h}};\mathbf{y})\prod\limits_{i=1}^{n-1}(x_{h}-y_{i})\prod\limits_{\begin{subarray}{c}i=1\\\
i\neq
h\end{subarray}}^{n-\lambda_{1}}(z_{2}-x_{i})\prod\limits_{i=3}^{(\lambda_{1}-\lambda_{2}+1)}\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq h\end{subarray}}^{n-\lambda_{1}+1}(z_{i}-x_{m})$
$\displaystyle=\mathfrak{S}_{\lambda^{\prime}}^{n-1}(\sigma^{\lambda_{1}-\lambda_{2}+1}(\mathbf{z});\mathbf{x}_{\hat{h}};\mathbf{y})\prod\limits_{i=1}^{n-1}(x_{h}-y_{i})\prod\limits_{i=2}^{(\lambda_{1}-\lambda_{2}+1)}\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq
h\end{subarray}}^{n-\lambda_{1}}(z_{i}-x_{m})\prod\limits_{i=3}^{\lambda_{1}-\lambda_{2}+1}(z_{i}-x_{n+1-\lambda_{1}})$
$\displaystyle=\mathfrak{S}_{(\lambda_{1},\lambda^{\prime})}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})|_{z_{1}=x_{h}}\prod\limits_{i=3}^{\lambda_{1}-\lambda_{2}+1}(z_{i}-x_{n+1-\lambda_{1}})=LHS|_{z_{1}=x_{h}}.$

∎

######  Proposition A.5.

Choose $\lambda\in\operatorname{Val}(n)$ with
$\operatorname{mul}(\lambda)=b>1$. Write
$\lambda=((\lambda_{1})^{b},\tilde{\lambda})$ for some $\tilde{\lambda}$. We
have

(A.8) $\displaystyle\mathfrak{S}_{((\lambda_{1})^{b},\tilde{\lambda})}^{n}$
$\displaystyle(\mathbf{z};\mathbf{x};\mathbf{y})\prod\limits_{i=1}^{b-1}(z_{i}-y_{n+1-b})\prod\limits_{i=b+2}^{b+\lambda_{1}-\tilde{\lambda}_{1}}(z_{i}-x_{n+1-\lambda_{1}})$
$\displaystyle=\partial^{\mathbf{z}}_{b}\left(\mathfrak{S}_{((\lambda_{1})^{b-1},\lambda_{1}-1,\tilde{\lambda})}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})\prod\limits_{i=1}^{\operatorname{mul}((\lambda_{1}-1,\tilde{\lambda}))-1}(z_{b+1}-y_{n+1-b-i})\right)$

###### Proof.

We view both sides as polynomials in $z_{1}$ and analyze their degrees to use
interpolation. By (5.1),
$\mathfrak{S}_{((\lambda_{1})^{b},\tilde{\lambda})}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})$
is of degree at most $(n-\lambda_{1}-b)$ in $z_{1}$ so the right hand side is
of degree at most $(n-\lambda_{1}-b+1)$ in $z_{1}$. Likewise,
$\mathfrak{S}_{((\lambda_{1})^{b-1},\lambda_{1}-1,\tilde{\lambda})}^{n}(\mathbf{z};\mathbf{x};\mathbf{y})$
is of degree at most $(n-\lambda_{1}-b+1)$ in $z_{1}$ so the right hand side
is of degree at most $(n-\lambda_{1}-b+1)$ in $z_{1}$. Since $b>1$, it is
enough to show that they coincide when we plug in $z_{1}=x_{a}$ for $1\leq
a\leq n-\lambda_{1}$.

We use induction on $b$. First assume $b=2$. Setting $z_{1}:=x_{a}$ on both
sides, A.3 gives

$\displaystyle LHS|_{z_{1}=x_{a}}$
$\displaystyle=\mathfrak{S}_{(\lambda_{1},\tilde{\lambda})}^{n-1}(\mathbf{z};\mathbf{x}_{\hat{a}};\mathbf{y})\prod\limits_{l=1}^{n-1}(x_{a}-y_{l})\prod\limits_{i=4}^{2+\lambda_{1}-\tilde{\lambda}_{1}}(z_{i}-x_{n+1-\lambda_{1}})$
$\displaystyle RHS|_{z_{1}=x_{a}}$
$\displaystyle=\partial^{\mathbf{z}}_{2}\left(\mathfrak{S}_{(\lambda_{1}-1,\tilde{\lambda})}^{n-1}(\sigma(\mathbf{z});\mathbf{x}_{\hat{a}};\mathbf{y})\prod\limits_{l=1}^{n-1}(x_{a}-y_{l})\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq
a\end{subarray}}^{n-\lambda_{1}}(z_{2}-x_{m})\prod\limits_{i=1}^{\operatorname{mul}((\lambda_{1}-1,\tilde{\lambda}))-1}(z_{3}-y_{n-1-i})\right).$

So the equality $LHS|_{z_{1}=x_{a}}=RHS|_{z_{1}=x_{a}}$ comes from A.4, with
the variables shifted by $\mathbf{x}\rightarrow\mathbf{x}_{\hat{a}}$ and
$\mathbf{z}\rightarrow\sigma(\mathbf{z})$.

For $b>2$, assume the statement holds for $(b-1)$. Setting $z_{1}:=x_{a}$ on
both sides, A.3 gives

$\displaystyle LHS|_{z_{1}=x_{a}}$
$\displaystyle=\mathfrak{S}_{((\lambda_{1})^{b-1},\tilde{\lambda})}^{n-1}(\mathbf{z};\mathbf{x}_{\hat{a}};\mathbf{y})\prod\limits_{l=1}^{n-b+1}(x_{a}-y_{l})\prod\limits_{i=2}^{b-1}(z_{i}-y_{n+1-b})\prod\limits_{i=b+2}^{b+\lambda_{1}-\tilde{\lambda}_{1}}(z_{i}-x_{n+1-\lambda_{1}})$
$\displaystyle RHS|_{z_{1}=x_{a}}$
$\displaystyle=\partial^{\mathbf{z}}_{b}\left(\mathfrak{S}_{((\lambda_{1})^{b-2},\lambda_{1}-1,\tilde{\lambda})}^{n-1}(\mathbf{z};\mathbf{x}_{\hat{a}};\mathbf{y})\prod\limits_{l=1}^{n-b+1}(x_{a}-y_{l})\prod\limits_{i=1}^{\operatorname{mul}((\lambda_{1}-1,\tilde{\lambda}))-1}(z_{b+1}-y_{n+1-b-i})\right).$

The equality $LHS|_{z_{1}=x_{a}}=RHS|_{z_{1}=x_{a}}$ comes from the induction
hypothesis with the variables shifted by
$\mathbf{x}\rightarrow\mathbf{x}_{\hat{a}}$ and
$\mathbf{z}\rightarrow\sigma(\mathbf{z})$. ∎

## References

  * [AGS20] Erik Aas, Darij Grinberg, and Travis Scrimshaw.  Multiline queues with spectral parameters.  Comm. Math. Phys., 374(3):1743–1786, 2020.
  * [AL14] Arvind Ayyer and Svante Linusson.  An inhomogeneous multispecies TASEP on a ring.  Adv. in Appl. Math., 57:21–43, 2014.
  * [AL18] Erik Aas and Svante Linusson.  Continuous multi-line queues and TASEP.  Ann. Inst. Henri Poincaré D, 5(1):127–152, 2018.
  * [AM13] Chikashi Arita and Kirone Mallick.  Matrix product solution of an inhomogeneous multi-species TASEP.  J. Phys. A, 46(8):085002, 11, 2013.
  * [AS16] Erik Aas and Jonas Sjöstrand.  A product formula for the TASEP on a ring.  Random Structures Algorithms, 48(2):247–259, 2016.
  * [BB93] Nantel Bergeron and Sara Billey.  RC-graphs and Schubert polynomials.  Experiment. Math., 2(4):257–269, 1993.
  * [Can16] Luigi Cantini.  Inhomogenous multispecies TASEP on a ring with spectral parameters, 2016\.  arXiv:1602.07921.
  * [Can17] Luigi Cantini.  Asymmetric simple exclusion process with open boundaries and Koornwinder polynomials.  Ann. Henri Poincaré, 18(4):1121–1151, 2017.
  * [CdGW15] Luigi Cantini, Jan de Gier, and Michael Wheeler.  Matrix product formula for Macdonald polynomials.  J. Phys. A, 48(38):384001, 25, 2015.
  * [CFF18] Rocco Chirivi’, Xin Fang, and Ghislain Fourier.  Degenerate Schubert varieties in type A, 2018.  arXiv:1808.01594.
  * [CMW18] Sylvie Corteel, Olya Mandelshtam, and Lauren Williams.  From multiline queues to Macdonald polynomials via the exclusion process, 2018.  arXiv:1811.01024.
  * [CW11] Sylvie Corteel and Lauren K. Williams.  Tableaux combinatorics for the asymmetric exclusion process and Askey-Wilson polynomials.  Duke Math. J., 159(3):385–415, 2011.
  * [CW18] Sylvie Corteel and Lauren K. Williams.  Macdonald-Koornwinder moments and the two-species exclusion process.  Selecta Math. (N.S.), 24(3):2275–2317, 2018.
  * [FI] OEIS Foundation Inc. (2019).  The On-Line Encyclopedia of Integer Sequences.  http://oeis.org.
  * [FK96] Sergey Fomin and Anatol N. Kirillov.  The Yang-Baxter equation, symmetric functions, and Schubert polynomials.  In Proceedings of the 5th Conference on Formal Power Series and Algebraic Combinatorics (Florence, 1993), volume 153, pages 123–143, 1996.
  * [FM07] Pablo A. Ferrari and James B. Martin.  Stationary distributions of multi-type totally asymmetric exclusion processes.  Ann. Probab., 35(3):807–832, 2007.
  * [KM05] Allen Knutson and Ezra Miller.  Gröbner geometry of Schubert polynomials.  Ann. of Math. (2), 161(3):1245–1318, 2005.
  * [LW12] Thomas Lam and Lauren Williams.  A Markov chain on the symmetric group that is Schubert positive?  Exp. Math., 21(2):189–192, 2012.
  * [Man01] Laurent Manivel.  Symmetric functions, Schubert polynomials and degeneracy loci, volume 6 of SMF/AMS Texts and Monographs.  American Mathematical Society, Providence, RI; Société Mathématique de France, Paris, 2001.  Translated from the 1998 French original by John R. Swallow, Cours Spécialisés [Specialized Courses], 3.
  * [Man17] Olya Mandelshtam.  Toric tableaux and the inhomogeneous two-species TASEP on a ring.  2017\.  arXiv:1707.02663.
  * [USW04] Masaru Uchiyama, Tomohiro Sasamoto, and Miki Wadati.  Asymmetric simple exclusion process with open boundaries and Askey-Wilson polynomials.  J. Phys. A, 37(18):4985–5002, 2004.


# Synthesizing impurity clustering in the edge plasma of tokamaks using neural
networks

Z. Lin  Aix-Marseille Université, CNRS, I2M, UMR 7373, 13453 Marseille, France
T. Maurel–Oujia  Aix-Marseille Université, CNRS, I2M, UMR 7373, 13453
Marseille, France  B. Kadoch  Aix-Marseille Université, CNRS, IUSTI, UMR 7343,
13453 Marseille, France  P. Krah  Aix-Marseille Université, CNRS, I2M, UMR
7373, 13453 Marseille, France  N. Saura  Aix-Marseille Université, CNRS, PIIM,
UMR 7345, 13397 Marseille, France  S. Benkadda  Aix-Marseille Université,
CNRS, PIIM, UMR 7345, 13397 Marseille, France  K. Schneider
<EMAIL_ADDRESS>Aix-Marseille Université, CNRS, I2M, UMR 7373, 13453
Marseille, France

###### Abstract

This work investigates the behavior of impurities in edge plasma of tokamaks
using high-resolution numerical simulations based on Hasegawa–Wakatani
equations. Specifically, it focuses on the behavior of inertial particles,
which has not been extensively studied in the field of plasma physics. Our
simulations utilize one-way coupling of a large number of inertial point
particles, which model plasma impurities. We observe that with Stokes number
($St$) which characterizes the inertia of particles being much less than one,
such light impurities closely track the fluid flow without pronounced
clustering. For intermediate $St$ values, distinct clustering appears, with
larger Stokes values, i.e. heavy impurities even generating more substantial
clusters. When $St$ is significantly large, very heavy impurities tend to
detach from the flow and maintain their trajectory, resulting in fewer
observable clusters and corresponding to random motion. A core component of
this work involves machine learning techniques. Applying three different
neural networks - Autoencoder, U-Net, and Generative Adversarial Network (GAN)
- to synthesize preferential concentration fields of impurities, we use
vorticity as input and predict impurity number density fields. GAN outperforms
the two others by aligning closely with direct numerical simulation data in
terms of probability density functions of the particle distribution and energy
spectra. This machine learning technique holds the potential to reduce
computational costs by eliminating the need to track millions of particles
modeling impurities in simulations.

††preprint: AIP/123-QED

##  I Introduction

The pursuit of a sustainable and abundant energy source has sparked greater
attention towards nuclear fusion, which is the process responsible for
powering celestial bodies such as stars. Nuclear fusion is distinct from
nuclear fission as it involves the fusion of light atomic nuclei to create
heavier ones, resulting in the release of a significant amount of energy.
Fusion is highly regarded as the ultimate goal in energy production due to its
capacity to generate significant power while causing minimal harm to the
environment.

Fusion reactors confine and heat a plasma, consisting of ions and electrons,
to enable fusion reactions at high temperatures. However, the task of
confining plasma poses a considerable challenge due to the presence of
multiple instabilities that can result in energy loss. Impurity accumulation
in the plasma core can result in heat loss through radiation, leading to a
decrease in confinement quality. Therefore, it is crucial to conduct a
thorough investigation of impurity in fusion plasma.

The dynamics of magnetically confined plasma flow are influenced by drift-wave
turbulence and zonal flows. The numerical high-resolution simulations for
modeling the plasma flow carried out within this work are based on the
Hasegawa–Wakatani model, which governs cross-field transport by electrostatic
drift waves in magnetically confined plasmas. A modified version of the
modelPushkarev, Bos, and Nazarenko (2013) is likewise investigated to take
into account zonal flows. The electric field perpendicular to magnetic field
lines is particularly significant because it strongly drives cross-field
fluxes, impacting edge pressure profiles and stability.Wang  _et al._ (2020);
Zhang, Guo, and Diamond (2020) $\bm{E}\times\bm{B}$ flows strongly influence
the motion of coherent structures, which can account for substantial particle
losses in tokamaks.Carralero  _et al._ (2017); D’Ippolito, Myra, and Zweben
(2011); Kube  _et al._ (2018) The resulting fluxes on walls can cause
sputtering, erosion, and impurity injection, further degrading confinement and
safe operation.D’Ippolito, Myra, and Zweben (2011); Kuang  _et al._ (2018,
2020)

For the particles, there has been growing interest in Lagrangian perspective
in recent years.Gheorghiu, Militello, and Rasmussen (2024); Kadoch  _et al._
(2022); Bos  _et al._ (2010) This approach involves the analysis of
trajectories of numerous tracers. Numerical simulations involve solving
equations to determine the trajectories of test particles within a specific
velocity field, such as the $\bm{E}\times\bm{B}$ velocity field. The
Lagrangian approach highlights the significant impact of coherent structures
on transport. The combined effect of eddy trapping and zonal shear flows often
leads to non-diffusive transport. Krasheninnikov, D’Ippolito, and Myra (2008);
del Castillo-Negrete, Carreras, and Lynch (2004); van Milligen, Sanchez, and
Carreras (2004) In fusion plasma research, previous studies have primarily
focused on passive flow tracers without considering inertial effects in the
edge plasma using the Hasegawa–Wakatani model.Futatani  _et al._ (2008a);
Futatani, Benkadda, and del Castillo-Negrete (2009); Futatani  _et al._
(2008b) We increase the complexity of these models and add the effect of
inertia, which has not been done so far, to the best of our knowledge.
Therewith we can study the behavior of heavy and light atoms and see the
impact on the impurity distribution. A crucial aspect of understanding
inertial impurities is their self-organization, such as clustering and void
formation, which can be quantified mathematically by deviation from Poissonian
statistics.Oujia, Matsuda, and Schneider (2020) Their dynamics can then be
analyzed using tesselation based techniques.Maurel-Oujia, Matsuda, and
Schneider (2024)

There is a growing interest in studying spatial patterns in turbulence
exploring machine learning.Ishikawa  _et al._ (2022); Jajima  _et al._ (2023)
Our study uses machine learning to estimate mesoscale inertial particle
clustering in turbulence, relying on flow field data rather than individual
particle information. We apply machine learning tools that have been
successfully implemented in the field of hydrodynamics, see Ref. Oujia  _et
al._ , 2022 and Maurel-Oujia  _et al._ , 2023. Three neural networks —
Autoencoder, U-Net, and Generative Adversarial Network (GAN) — are trained to
predict impurity particle densities from vorticity field snapshots. The
results are visually and statistically compared with actual DNS particle data.
This approach reduces the cost for predicting particle density distributions.

The paper is structured as follows: Sec. II outlines the theoretical basis
used for simulation, introducing the Hasegawa–Wakatani edge plasma model and
the motion of inertial impurity particles. Sec. III details the numerical
simulations and results. In Sec. IV, we explore machine learning models such
as Autoencoders, U-net, and GAN to synthesize impurity concentration fields.
Finally, Sec. V summarizes our findings and suggests potential directions for
future investigations.

##  II Models for Edge Plasma Turbulence and impurity particles

###  II.1 Hasegawa–Wakatani model

The numerical simulations are based on Hasegawa–Wakatani (HW) equations for
plasma edge turbulence driven by drift-wave instability.Hasegawa and Wakatani
(1983) Our focus is on the two-dimensional, slab geometry of the HW model, see
e.g. Ref. Kadoch  _et al._ , 2022. Fig. 1 depicts a representation of the flow
configuration. The magnetic field lines are assumed to be constant, straight
and perpendicular to the slab. Ions are considered cold, and the effects of
the temperature gradient are ignored.

Figure 1: Illustration of the two-dimensional slab geometry in the tokamak
edge used in the Hasegawa–Wakatani system. The radial direction is represented
by $x$, while $y$ is the poloidal direction. There is an imposed mean plasma
density gradient $\nabla{n_{0}}$ in the radial direction. The 2D flow is
computed within a domain whose size is equivalent to 64 times the Larmor
radius ($\rho_{s}$). This computational domain is marked with a green square.
On the right, a vorticity field for the classical Hasegawa–Wakatani (cHW) is
shown. The figure is adapted from Ref. Bos  _et al._ , 2010.

.

The HW model consists of two partial differential equations that describe the
time evolution of the plasma electrostatic potential $\phi$ and fluctuating
plasma density $n$:

$\left(\frac{\partial}{\partial
t}-\nu\nabla^{2}\right)\nabla^{2}\phi=\left[\nabla^{2}\phi,\phi\right]+c(\phi-n),$
(1)

$\left(\frac{\partial}{\partial
t}-D\nabla^{2}\right)n=[n,\phi]-\kappa\frac{\partial\phi}{\partial
y}+c(\phi-n),$ (2)

All quantities are dimensionless and normalized as described in Refs. Bos  _et
al._ , 2010; Futatani  _et al._ , 2011. The constant $D$ and $\nu$ are cross-
field diffusion coefficient of plasma density fluctuations and kinematic
viscosity, respectively. The constant $\kappa$ defined as
$\kappa\equiv-\partial_{x}\ln n_{0}$, is a measure of the density gradient.
The Poisson bracket is defined as: $[A,B]=\frac{\partial A}{\partial
x}\frac{\partial B}{\partial y}-\frac{\partial A}{\partial y}\frac{\partial
B}{\partial x}$. In this system, the electrostatic potential $\phi$ acts as a
stream-function for the $\bm{E}\times\bm{B}$ flow velocity, denoted by
$\bm{u}=\nabla^{\perp}\phi$, where
$\nabla^{\perp}=\left(-\partial_{y},\partial_{x}\right)$. Thus, we have
$u_{x}=-\partial\phi/\partial y$ and $u_{y}=\partial\phi/\partial x$. The
vorticity is defined as $\omega=\nabla^{2}\phi$. The turnover time of
turbulent eddies, $\tau_{\eta}$, is defined as $1/\sqrt{2Z_{ms}}$, with
$Z_{ms}$ denoting the mean-square vorticity. The adiabaticity parameter $c$
is:

$c=\frac{T_{e}k_{\parallel}^{2}}{\mathrm{e}^{2}n_{0}\eta\omega_{ci}}$ (3)

where $\eta$ is the electron resistivity. Here $k_{\parallel}$ is the
effective parallel wavenumber. The parameter $c$ controls the phase difference
between the electrostatic potential and the plasma density fluctuations. The
case $c\gg 1$ is known as the adiabatic limit. In this limit, the
Hasegawa–Wakatani model reduces to the Hasegawa–Mima equation, and the
electrons in the plasma follow a Boltzmann distribution. In that case, there
is no phase shift between $\phi$ and $n$. Conversely, when $c\ll 1$ (the
hydrodynamic limit), the system reduces to a form that is analogous to a two-
dimensional Navier–Stokes equation. In this regime, the density fluctuations
are passively advected by the $\bm{E}\times\bm{B}$ flow.Kadoch  _et al._
(2022) The most physically interesting situation in fusion plasma is when
$c\sim 1$, known as the quasi-adiabatic regime. In this regime, there is a
phase shift between $\phi$ and $n$. This phase shift permits the presence of
particle transport.

In the context of the classical Hasegawa–Wakatani model (cHW), as detailed
above, zonal flows are not observed. To specifically obtain zonal flows, a
modified Hasegawa–Wakatani model (mHW) can be considered. This modification
was introduced in Ref. Pushkarev, Bos, and Nazarenko, 2013 where the coupling
term $c(\phi-n)$ is set to zero for modes $k_{y}=0$.

###  II.2 Impurity particles model

Previous studies Futatani  _et al._ (2008a); Futatani, Benkadda, and del
Castillo-Negrete (2009); Futatani  _et al._ (2008b) have assumed that impurity
particles perfectly follow the fluid flow without deviating from the flow
path. However, this assumption may not always be valid, particularly when the
particles have significant mass, resulting in observable inertial effects.

####  II.2.1 Stokes number

To account for the inertial effects, the concept of the Stokes number ($St$)
is introduced. Maxey (1987) The Stokes number is a dimensionless parameter
used to quantify the inertia of particles in a fluid flow. The Stokes number
is defined as:

$St=\frac{\tau_{p}}{\tau_{\eta}},$ (4)

where $\tau_{{p}}$ is the impurity particle’s relaxation time. This is a
measure of how quickly an impurity particle can respond to changes in the
fluid flow. It is determined by factors like the particle’s size, mass, as
well as the viscosity of the fluid. The turbulent eddy turnover time is
$\tau_{\eta}$. This is a measure of how quickly the fluid flow is changing. A
high Stokes number suggests that the particle’s motion is dominated by its
inertia, meaning it will continue moving in its current direction even when
the fluid flow changes. Note that the Stokes number is a simplification, and
real-world particle-fluid interactions may involve other complex phenomena,
like particle-particle interactions, two-way coupling, finite particle-size
effects, etc. Nevertheless, the Stokes number provides a valuable prediction
of particle behavior in a turbulent flow.

####  II.2.2 Motion of inertial impurity particles

We assume the behavior of impurities does not influence the plasma’s overall
dynamics. This study is focused on tracking individual particles using
Lagrangian mechanics. The equations governing the motion of these particles
are:

$\frac{d\bm{x}_{imp,j}}{dt}=\bm{v}_{imp,j}$ (5)

$\frac{d\bm{v}_{imp,j}}{dt}=-\frac{\bm{v}_{imp,j}-\bm{u}_{imp,j}}{\tau_{p}}$
(6)

The particles are denoted by $j$ ($j=1,\ldots,N$, where $N$ is the total
number of impurity particles in the system). Here, $\bm{v}_{imp,j}$ is the
velocity vector of the $j$-th impurity particle, and $\bm{u}_{imp,j}$
represents the fluid velocity vector at the location of the particle,
$\bm{x}_{imp,j}$. The force resisting the motion of these particles is assumed
to be in the form of Stokes drag, which is directly proportional to the
relative velocity between the particle and the fluid. The constant of
proportionality, $1/\tau_{p}$, is the inverse of particle relaxation
time.Onishi, Baba, and Takahashi (2011); Obligado  _et al._ (2014); Matsuda
and Onishi (2019) In the initial phase of our study, we excluded
electromagnetic forces to focus on establishing a fundamental understanding
without introducing additional complexity. The effects of electromagnetic
forces on the particles will be integrated incrementally in future studies.

##  III Numerical Simulations

###  III.1 Flow configurations and physical parameters

The simulation was conducted within a domain with periodic boundary condition
that spans an area of $64^{2}$ and it was discretized using a resolution of
$1024^{2}$ grid points. The time step was $5\times 10^{-4}$, while values for
diffusivity of plasma density $D$ and kinematic viscosity $\nu$ were both set
to $5\times 10^{-3}$. The value for $\kappa$ was set to 1. Simulations for the
flow begin with Gaussian random initial conditions and run until achieving a
saturated, fully developed turbulent flow (see Ref. Kadoch  _et al._ , 2022).
After this phase,$10^{6}$ uniformly distributed inertial impurity particles
are injected.

###  III.2 Numerical schemes

The systems of HW equations, specifically Eq. (1) and Eq. (2) are solved using
a Fourier pseudospectral method. It is particularly well-suited for problems
with periodic boundary conditions and provides excellent accuracy for smooth
solutions.Canuto  _et al._ (2007) The equations governing the motion of the
impurity particles, namely Eq. (5) and Eq. (6), are solved using a second-
order Runge–Kutta (RK2) scheme. This method provides a balance between
accuracy and computational efficiency. Simulating a large ensemble of
particles presents a considerable computational challenge. Our solution
leverages High-Performance Computing (HPC) to overcome this complexity,
employing parallel computing with Message Passing Interface (MPI). For further
details on the numerical method and its implementation, we refer to Ref.
Kadoch  _et al._ , 2022.

###  III.3 Results

Figure 2: Vorticity fields without impurity particles in fully developed
turbulence regime.

Figure 3:  Kinetic energy spectra of the flow. The slope $k^{-3.5}$ is
provided for reference.

Fig. 2 shows the vorticity field (without impurity particles) in a fully
developed turbulence regime of classical Hasegawa–Wakatani model (cHW) for
different $c$ values and modified Hasegawa–Wakatani model (mHW) for $c=2$.
Larger vortices were exhibited in hydrodynamic cases (cHW, $c=0.01$), while
smaller ones were observed in quasi-adiabatic (cHW, $c=0.7$) and adiabatic
cases (cHW, $c=2$). For the modified Hasegawa–Wakatani model, with $c=2$,
zonal flows are presented. To quantify the distribution of kinetic energy of
the flow across different scales of the turbulent flow, we compute the Fourier
energy spectra:Bassenne  _et al._ (2017)

$E(k)=\frac{1}{2}\,\sum_{k-\frac{1}{2}\leq|{\bm{k}}|<k+\frac{1}{2}}|\widehat{\bm{u}}({\bm{k}})|^{2}$
(7)

where
$\widehat{\bm{u}}({\bm{k}})=\int\int\,{\bm{u}}({\bm{x}})\exp(-i2\pi{\bm{k}}\cdot{\bm{x}})\,d{\bm{x}}$
is the 2D Fourier transform of the flow velocity, $i=\sqrt{-1}$ and
${\bm{k}}=(k_{x},k_{y})$ the wave vector. The summation in Eq. (7) is carried
out over concentric shells in wave number space. In Fig. 3 the Fourier energy
spectra $E(k)$ are shown. We observe a peak at smaller wavenumbers for the
cases $c=0.01$ (cHW and mHW) and $c=2$ (mHW), indicating the presence of
larger structures. Furthermore, the spectra display a power-law scaling, with
an exponent nearing -3.5 for the case $c=0.7$.

In this paper, we focus on the quasi-adiabatic regime ($c=0.7$), which is the
relevant case for edge plasma of tokamaks. For readers interested in zonal
flow, relevant discussions can be found in Appendix A. Fig. 4 shows the
visualization of vorticity fields and $10^{4}$ superimposed impurity particles
for various Stokes numbers in quasi-adiabatic regime. From Fig. 4 we can
observe clustering of impurity particles for the case of $St=0.25$ and $St=1$,
while for $St=0$ and $St=50$, the impurity density field homogeneously fills
the entire physical domain and does not show significant spatial correlation.
When the particles are very light, the response time of the particles is much
smaller than the turnover time of the turbulent eddies, i.e.,
$\tau_{p}\ll\tau_{\eta}$ ($St\sim 0$). Their response time is so quick that
they can adjust their velocity to match any changes in the flow immediately.
They behave as usual tracers like those in previous studies. Futatani  _et
al._ (2008a); Futatani, Benkadda, and del Castillo-Negrete (2009); Futatani
_et al._ (2008b) When $\tau_{p}\gg\tau_{\eta}$ ($St\gg 1$), the inertia of the
impurity particles is so high that they are almost unaffected by the eddies.
Particles are moving ballistically and are randomly distributed, resulting in
less clustering and more dispersion. When the response time of the particles
is comparable to the eddy turnover time, i.e., $\tau_{p}\sim\tau_{\eta}$
($St\sim 1$), particles are subject to centrifugal effect of turbulent eddies.
Eddies can accelerate the particles to a degree where they are thrown out of
the eddy due to their inertia. This ejection and the subsequent movement
towards areas of lower vorticity result in particle clustering. The clustering
of impurity particles results in inhomogeneous poloidal distributions ($y$
direction in the domain) of the impurity density, which can modify radial
collisional transport.Angioni (2021)

Figure 4: Visualization of vorticity fields in fully developed turbulence
regime and superimposed impurity particles ($10^{4}$ impurity particles) for
various Stokes numbers in the case of $c=0.7$ (cHW).

To analyze and statistically interpret the distribution of impurity particles,
we increase the number of particles from $10^{4}$ to $10^{6}$. Displaying the
position of each particle in the flow is impractical due to this large
quantity. Instead, we visualize the particle density. We achieve this by
dividing the entire domain into $512\times 512$ boxes and defining the
impurity density $n_{imp}$ as the number of particles in each box, i.e. we
compute a histogram. Fig. 5 shows the impurity density field for the case of
$c=0.7$ considering different Stokes numbers. When $St=0$, the impurity
density is distributed uniformly throughout the flow. As $St$ increases,
variations in impurity density become apparent, with some areas exhibiting
high density and others low. This uneven distribution indicates the formation
of voids and clustering of particles. When $St=50$, the impurity density
returns to a more uniform distribution. The impurity density fields exhibit
patterns similar to those seen in the earlier visualization involving $10^{4}$
impurity particles.

To analyze the impurity density distribution at different scales, we compute
the Fourier energy spectra similar to Eq. (7) replacing the flow velocity
${\bm{u}}$ by the impurity density distribution $n_{imp}$. The spectral
characteristics of impurity density field are shown in Fig. 6 in terms of the
Fourier energy spectra. The variance of the impurity density is equivalent to
the cumulative sum of the corresponding spectral energy. This is visually
represented as the area under the energy spectra curve. As observed, the
variance initially increases with the Stokes number up to $St=1$, indicating
stronger fluctuations. From $St=1$ to $St=50$, there is a decrease, suggesting
less fluctuations. For $St=0$, representing tracer particles without inertia
as in previous studies, the impurity particles are distributed randomly.
Random spatial distribution of particles follows a Poisson probability density
function (PDF) with mean and variance equal to the mean density $\langle
n_{imp}\rangle$, see e.g. Ref. Bassenne  _et al._ , 2017. The slope of the
spectra for $St=0$ in 2D is 1, as prescribed by the equidistribution of
variance among all wavenumbers. The spectrum for $St=50$ is observed to
collapse on those of $St=0$ at small scales (high wavenumbers), thus
resembling the same characteristics as for randomly distributed particles.

Figure 5: Visualization of impurity density fields ($10^{6}$ impurity
particles) for various Stokes numbers in the classical Hasegawa–Wakatani model
($c=0.7$).

Figure 6: Energy spectra of impurity density for $St$ = 0, 0.25, 1 and 50 in
the case of $c=0.7$(cHW). The slope $k$ is provided for reference.

##  IV Neural Networks for Synthesizing Preferential Concentration of
Particles

Simulating flow without impurity particles to obtain the vorticity field is
not computationally expensive, but including and tracking $10^{6}$ impurity
particles increases costs. Since a specific statistically stationary
distribution of impurity density is linked to a specific vorticity field, we
propose creating a neural network that can directly predict the distribution
of impurity density based on the vorticity field. This network uses as input
the vorticity field (easily obtained from DNS) and produces as output the
corresponding impurity density field for a specific Stokes number. With this
neural network, we can run DNS to get the vorticity field, use it as input for
the network, and efficiently predict the impurity density distribution for a
specific Stokes number. This approach eliminates the need for costly
simulations involving millions of particles.

We will evaluate and compare three different neural network architectures that
have previously been utilized in particle-laden hydrodynamic turbulence
studies:Oujia  _et al._ (2022); Maurel-Oujia  _et al._ (2023) the Autoencoder,
U-Net, and Generative Adversarial Network (GAN). To determine the most
effective neural network architecture for accurately predicting the impurity
density fields, a thorough analysis of the statistical attributes of the
impurity density fields generated by these neural networks is carried out. In
the following, we will consider the cHW model for the quasi-adiabatic case
($c=0.7$).

###  IV.1 Different machine learning models: Autoencoder, U-Net, and GAN

Figure 7: The neural network architectures are illustrated as: Autoencoder in
purple; U-Net in purple and red; GAN in purple, red, and green.

The models used in this study are based on those utilized in Ref. Oujia  _et
al._ , 2022 and Maurel-Oujia  _et al._ , 2023. Fig. 7 displays the
architectures tailored from the standard structure, which we have implemented
using TensorFlow.

  * •

Autoencoder: Autoencoder Rumelhart, Hinton, and Williams (1986) is the basic
model in our study, comprising of an encoder that compresses the data, and a
decoder that reconstructs the output from this compressed version.

  * •

U-Net: U-Net Ronneberger, Fischer, and Brox (2015) is essentially an
Autoencoder but with added skip connections. These connections facilitate non-
sequential connections between layers which helps in the preservation of
information at different scales throughout the network.

  * •

GAN: GAN (Generative Adversarial Network) consists of two neural
networks,Goodfellow  _et al._ (2020) a Generator and a Discriminator, trained
simultaneously through adversarial processes. The Generator attempts to
produce synthetic data, while the Discriminator tries to distinguish between
real and synthetic data. We use U-Net as the generator.Isola  _et al._ (2017)
More details can be found in Ref. Oujia  _et al._ , 2022 and Maurel-Oujia  _et
al._ , 2023.

In this study, we use the data of one long simulation in the statistically
stationary regime. This simulation has a temporal span equivalent to 220 eddy
turnover times. Within this simulation, we have captured 400 snapshots of the
vorticity field. For each snapshot of the vorticity field, there are 5
corresponding impurity density fields ($St=0.01,0.05,0.25,1,5$). The neural
network is then trained to learn the relation between the vorticity field and
the corresponding impurity density field. 70% of the dataset is for training
the model and 30% for evaluating its performance. Vorticity data are
normalized between -1 and 1. Impurity density data are normalized between 0
and 1 to enable the use of binary cross-entropy.

The Adam optimizer Kingma and Ba (2014) was used to train the neural networks.
Training on a single sample from the training set is termed a ‘step’. The
Autoencoder and U-net models were trained for 28,000 training steps. The GAN
model required more steps for convergence and was trained for 56,000 steps The
Python codes were run on an NVIDIA graphics card Tesla V100 32GB.

###  IV.2 Visual Results

Fig. 8 shows the vorticity field and corresponding impurity density field
($St=1$) obtained from DNS for the case of $c=0.7$ (cHW). Fig. 9 illustrates
impurity density field ($St=1$) predicted from Autoencoder, Unet and GAN for
$c=0.7$ (cHW) using vorticity field as input. When comparing DNS impurity
density fields with those predicted by the Autoencoder and U-Net, we find that
both DNS and predicted fields exhibit similar structures, but the clustering
in the predicted fields is slightly fuzzy. The GAN model successfully predicts
void areas and filament-like structures to some degree of success, but it
sometimes struggles with predicting the details of regions with impurities.

Figure 8: Vorticity field (left)) and corresponding impurity density field
(right) ($St=1$) obtained from DNS for the case of $c=0.7$, cHW.

Figure 9: From left to right: For the case of $c=0.7$ (cHW), impurity density
field ($St=1$) predicted from Autoencoder, U-Net and GAN, respectively.

###  IV.3 Statistical analysis

To better understand the performance of our models, we analyze the results
using the probability density functions (PDFs) and energy spectra of the
impurity number density. Fig. 10 shows the PDFs of DNS data and predicted data
using three different architectures for $St=1$. From Fig. 10 we observe that
the GAN model excels in predicting the density distribution, as its
predictions closely align with DNS. In contrast, U-Net and Autoencoder exhibit
shortcomings, particularly in predicting higher densities. The deviations in
the PDFs of U-Net and the Autoencoder suggest that these models may struggle
to accurately replicate the intricate characteristics present in high-density
regions.

Fig. 11 shows the energy spectra of DNS data and predicted data using three
different architectures for $St=1$. In Fig. 11 we observe that the GAN’s
energy spectra are nearly identical to that of the DNS data. This demonstrates
the remarkable accuracy of the GAN model in predicting the energy spectra of
the impurity density. On the contrary, both U-Net and the Autoencoder exhibit
a significant drop at high wave numbers, which correspond to fine structures
in the data. This indicates that U-Net and the Autoencoder are unable to
accurately capture these fine-scale structures and may struggle to represent
the intricate details present in the impurity density.

Figure 10: PDFs of impurity density of DNS data and predicted data using three
different architectures for $St=1$.

Figure 11: Energy spectra of impurity density of DNS data and predicted data
using three different architectures for $St=1$.

Figure 12: PDFs of impurity density of DNS data (solid lines) and predicted
data using GAN (dashed lines) for different $St$.

Figure 13: Energy spectra of impurity density of DNS data (solid lines) and
predicted data using GAN (dashed lines) for different $St$.

Overall, the analysis of the PDFs and energy spectra reaffirms that GAN
outperforms the other two for $St=1$. Thus, we utilize GAN in the following to
investigate more Stokes numbers values. From the PDF of impurity density as
shown in Fig. 12, we can see that for $St=0.25$, GAN predicts lower impurity
density when compared with DNS. For small $St$ ($St=0.01,0.05$), GAN predicts
more particles than DNS. When dealing with a large Stokes number ($St=5$), the
PDF is similar to the DNS results. In the case of $St=1$, the GAN predictions
align perfectly with DNS. For energy spectra, as shown in Fig. 13, it is
noticeable that the GAN predictions align very well with the DNS. In the GAN
predicted energy spectra, notable spikes were identified at frequencies
corresponding to the edges and corners of the domain. This could be because
the initial encoding of domain periodicity was inadequate. To correct this,
adjustments in the code are needed for true periodic convolution.
Unfortunately, as of the time of this writing, TensorFlow does not yet provide
native support for periodic padding.

In our study, we train the AI model specifically for $c=0.7$, which is a
typical value for edge plasma. This model predicts impurity density fields for
$St=0.01,0.05,0.25,1$ and $5$ using the vorticity field at $c=0.7$ as input.
The model’s validity is specific to $c=0.7$ and the aforementioned Stokes
numbers. Its extrapolability is shown and discussed in Appendix B.

Now we know that GAN is the optimal model. Next, we can conduct simulations to
generate more vorticity fields without particles, a process which is
computationally efficient. We then input these vorticity fields into GAN to
predict the impurity density fields. Subsequently, we can perform detailed
analyses of transport characteristics and clustering, using methods like
multiresolution analysis and Voronoi tessellation. There are several possible
ways to enhance the model. One method is to adjust architectural features like
the number of layers, filter sizes, and channel counts. Another effective
approach is to integrate physical constraints into the model, ensuring that
predictions align with established physical laws, such as mass or energy
conservation, for more accurate and plausible results.

##  V Conclusions

This study focuses on inertial impurity within fusion plasmas, which is
different from previous literature where the focus has typically been on
‘passive tracers’ following the flow. We used high-resolution numerical
simulations based on Hasegawa–Wakatani equations with impurity particles
having different inertia. We studied the quasi-adiabatic regime ($c=0.7$),
typical for Tokamak edge plasma. Our findings show that the inertia of
particles, quantified by the Stokes number $St$, plays a vital role in their
behavior. At lower $St$ values, impurities move along the flow streamlines. As
$St$ increases, impurities tend to cluster at low-vorticity regions. At very
high $St$ values, impurities show random movement due to substantial inertia,
resulting in much less clustering.

As simulating $10^{6}$ or more impurity particles is a computationally
intensive task, we utilized machine learning techniques to create a surrogate
model. Three neural networks - Autoencoder, U-Net, and GAN - were developed to
synthesize impurity densities, using vorticity fields as input. All models
produced visually comparable results to simulations, with varying degrees of
success. The GAN closely matches DNS values in terms of probability density
functions and energy spectra of density. This confirms GAN’s effectiveness in
modeling both the general distribution and finer structures. Based on the
results, future work will involve: Improving GAN by adjusting the
architectures or involving physical constraints directly into the model;
running simulations for more vorticity fields, then using the GAN to predict
impurity fields for analyzing clustering and transport properties; Including
electromagnetic interactions for more accurate impurity particle models.

###### Acknowledgements.

ZL, TMO and KS acknowledge the financial support from I2M and the French
Federation for Magnetic Fusion Studies (FR-FCM) and the Eurofusion consortium,
funded by the Euratom Research and Training Programme under Grant Agreement
No. 633053. The views and opinions expressed herein do not necessarily reflect
those of the European Commission. ZL, TMO, BK, PK and KS acknowledge partial
funding from the Agence Nationale de la Recherche (ANR), project CM2E, grant
ANR-20-CE46-0010-01. Centre de Calcul Intensif d’Aix-Marseille is acknowledged
for providing access to its high performance computing resources.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Impurities in zonal flows

Figure 14: Visualization of vorticity fields in fully developed turbulence
regime and $10^{4}$ superimposed impurity particles for various Stokes numbers
in the modified Hasegawa–Wakatani model ($c=2$).

Figure 15: Visualization of impurity density fields ($10^{6}$ impurity
particles) for various Stokes numbers in the modified Hasegawa–Wakatani model
($c=2$).

In the modified Hasegawa–Wakatani model at $c=2$ with $10^{4}$ impurity
particles as shown in Fig. 14, the zonal flow is very strong and the vorticity
values are weaker compared to cHW. The particles appear to be distributed
evenly across the fluid, and in proximity to the zonal flow, it seems that
clustering is not evident To gain deeper insights, we examined the impurity
density fields obtained from $10^{6}$ impurity particles. As depicted in Fig.
15, there is some preferential concentration of impurities in the shear layer
region for $St=0.25$, 1 and 5. The impurity particles are ejected by the zonal
flow, and for St=1 the effect seems to be the strongest. For $St=0$ and 50 the
impurity particles are randomly distributed.

##  Appendix B Extrapolability

In the following we test the extrapolability of the GAN model for different
$c$-values. The model has been trained using the data with $c=0.7$. As input
data we then use vorticity fields for $c$-values different from 0.7 and assess
the impurity densities obtained with the GAN model. The results are shown from
Fig. 16 to Fig. 19. In the hydrodynamic regime ($c=0.01$), we observe that GAN
predicts the PDF and energy spectra of impurity density very well at $St=1$.
At other Stokes numbers, discrepancies arise between GAN predictions and DNS
results in terms of PDFs and energy spectra. In the adiabatic regime ($c=2$),
the PDFs and energy spectra of impurity density predicted by GAN tend to have
larger values than the DNS results. Note that the results are good from
statistical point of view, but the obtained impurity density fields differ
visually from the DNS results.

Figure 16: PDFs of impurity density of DNS data (solid lines) and predicted
data using GAN (dashed lines) for different $St$ in the hydrodynamic regime
($c=0.01$).

Figure 17: Energy spectra of impurity density of DNS data (solid lines) and
predicted data using GAN (dashed lines) for different $St$ in hydrodynamic
regime ($c=0.01$).

Figure 18: PDFs of impurity density of DNS data (solid lines) and predicted
data using GAN (dashed lines) for different $St$ in the adiabatic regime
($c=2$).

Figure 19: Energy spectra of impurity density of DNS data (solid lines) and
predicted data using GAN (dashed lines) for different $St$ in the adiabatic
regime ($c=2$).

## References

  * Pushkarev, Bos, and Nazarenko (2013) A. V. Pushkarev, W. J. T. Bos, and S. V. Nazarenko, “Zonal flow generation and its feedback on turbulence production in drift wave turbulence,” Phys. Plasmas 20, 042304 (2013).
  * Wang  _et al._ (2020) H. Wang, H. Guo, G. Xu, A. Leonard, X. Wu, M. Groth, A. Jaervinen, J. Watkins, T. Osborne, D. Thomas, D. Eldon, P. Stangeby, F. Turco, J. Xu, L. Wang, Y. Wang, and J. Liu, “First evidence of local $\mathbf{E}\times{}\mathbf{B}$ drift in the divertor influencing the structure and stability of confined plasma near the edge of fusion devices,” Phys. Rev. Lett. 124, 195002 (2020).
  * Zhang, Guo, and Diamond (2020) Y. Zhang, Z. Guo, and P. Diamond, “Curvature of radial electric field aggravates edge magnetohydrodynamics mode in toroidally confined plasmas,” Phys. Rev. Lett. 125, 255003 (2020).
  * Carralero  _et al._ (2017) D. Carralero, M. Siccinio, M. Komm, S. Artene, F. D’Isa, J. Adamek, L. Aho-Mantila, G. Birkenmeier, M. Brix, G. Fuchert, M. Groth, T. Lunt, P. Manz, J. Madsen, S. Marsen, H. Müller, U. Stroth, H. Sun, N. Vianello, M. Wischmeier, E. Wolfrum, A. Team, C. Team, J. Contributors, and T. Team, “Recent progress towards a quantitative description of filamentary SOL transport,” Nucl. Fusion 57, 056044 (2017).
  * D’Ippolito, Myra, and Zweben (2011) D. D’Ippolito, J. Myra, and S. Zweben, “Convective transport by intermittent blob-filaments: Comparison of theory and experiment,” Phys. Plasmas 18, 060501 (2011).
  * Kube  _et al._ (2018) R. Kube, O. Garcia, A. Theodorsen, D. Brunner, A. Kuang, B. LaBombard, and J. Terry, “Intermittent electron density and temperature fluctuations and associated fluxes in the alcator C-Mod scrape-off layer,” Plasma Phys. Control. Fusion 60, 065002 (2018).
  * Kuang  _et al._ (2018) A. Kuang, N. Cao, A. Creely, C. Dennett, J. Hecla, B. LaBombard, R. Tinguely, E. Tolman, H. Hoffman, M. Major, J. R. Ruiz, D. Brunner, P. Grover, C. Laughman, B. Sorbom, and D. Whyte, “Conceptual design study for heat exhaust management in the ARC fusion pilot plant,” Fusion Eng. Des. 137, 221–242 (2018).
  * Kuang  _et al._ (2020) A. Kuang, S. Ballinger, D. Brunner, J. Canik, A. Creely, T. Gray, M. Greenwald, J. Hughes, J. Irby, B. LaBombard, and et al, “Divertor heat flux challenge and mitigation in SPARC,” J. Plasma Phys. 86 (2020).
  * Gheorghiu, Militello, and Rasmussen (2024) T. Gheorghiu, F. Militello, and J. J. Rasmussen, “On the transport of tracer particles in two-dimensional plasma edge turbulence,” Phys. Plasmas 31 (2024).
  * Kadoch  _et al._ (2022) B. Kadoch, D. del Castillo-Negrete, W. J. T. Bos, and K. Schneider, “Lagrangian conditional statistics and flow topology in edge plasma turbulence,” Phys. Plasmas 29 (2022).
  * Bos  _et al._ (2010) W. J. T. Bos, B. Kadoch, S. Neffaa, and K. Schneider, “Lagrangian dynamics of drift-wave turbulence,” Physica D 239, 1269–1277 (2010).
  * Krasheninnikov, D’Ippolito, and Myra (2008) S. I. Krasheninnikov, D. A. D’Ippolito, and J. R. Myra, “Recent theoretical progress in understanding coherent structures in edge and SOL turbulence,” J. Plasma Phys. 74, 679–717 (2008).
  * del Castillo-Negrete, Carreras, and Lynch (2004) D. del Castillo-Negrete, B. A. Carreras, and V. E. Lynch, “Fractional diffusion in plasma turbulence,” Phys. Plasmas 11, 3854–3864 (2004).
  * van Milligen, Sanchez, and Carreras (2004) B. P. van Milligen, R. Sanchez, and B. A. Carreras, “Probabilistic finite-size transport models for fusion: Anomalous transport and scaling laws,” Phys. Plasmas 11, 2272–2285 (2004).
  * Futatani  _et al._ (2008a) S. Futatani, S. Benkadda, Y. Nakamura, and K. Kondo, “Multiscaling dynamics of impurity transport in drift-wave turbulence,” Phys. Rev. Lett. 100, 025005 (2008a).
  * Futatani, Benkadda, and del Castillo-Negrete (2009) S. Futatani, S. Benkadda, and D. del Castillo-Negrete, “Spatiotemporal multiscaling analysis of impurity transport in plasma turbulence using proper orthogonal decomposition,” Phys. Plasmas 16, 042506 (2009).
  * Futatani  _et al._ (2008b) S. Futatani, S. Benkadda, Y. Nakamura, K. Kondo, and S. Hamaguchi, “Anomalous scaling of impurity transport in drift wave turbulence,” Contrib. Plasma Phys. 48, 111–115 (2008b).
  * Oujia, Matsuda, and Schneider (2020) T. Oujia, K. Matsuda, and K. Schneider, “Divergence and convergence of inertial particles in high Reynolds number turbulence,” J. Fluid Mech. 905, A14 (2020).
  * Maurel-Oujia, Matsuda, and Schneider (2024) T. Maurel-Oujia, K. Matsuda, and K. Schneider, “Computing differential operators of the particle velocity in moving particle clouds using tessellations,” J. Comput. Phys. 498, 112658 (2024).
  * Ishikawa  _et al._ (2022) R. T. Ishikawa, M. Nakata, Y. Katsukawa, Y. Masada, and T. L. Riethmüller, “Multi-scale deep learning for estimating horizontal velocity fields on the solar surface,” Astron. Astrophys. 658, A142 (2022).
  * Jajima  _et al._ (2023) Y. Jajima, M. Sasaki, R. T. Ishikawa, M. Nakata, T. Kobayashi, Y. Kawachi, and H. Arakawa, “Estimation of 2d profile dynamics of electrostatic potential fluctuations using multi-scale deep learning,” Plasma Phys. Control. Fusion 65, 125003 (2023).
  * Oujia  _et al._ (2022) T. Oujia, S. S. Jain, K. Matsuda, K. Schneider, J. West, and K. Maeda, “Neural networks for synthesizing preferential concentration of particles in isotropic turbulence,” in  _Proceedings of the Summer Program 2022, Center for Turbulence Research_ (Stanford University, 2022) pp. 153–162.
  * Maurel-Oujia  _et al._ (2023) T. Maurel-Oujia, S. S. Jain, K. Matsuda, K. Schneider, J. West, and K. Maeda, “Neural network models for preferential concentration of particles in two-dimensional turbulence,” arXiv preprint arXiv:2312.14829 (2023).
  * Hasegawa and Wakatani (1983) A. Hasegawa and M. Wakatani, “Plasma edge turbulence,” Phys. Rev. Lett. 50, 682–686 (1983).
  * Futatani  _et al._ (2011) S. Futatani, W. J. T. Bos, D. del Castillo-Negrete, K. Schneider, S. Benkadda, and M. Farge, “Coherent vorticity extraction in resistive drift-wave turbulence: Comparison of orthogonal wavelets versus proper orthogonal decomposition,” C. R. Phys. 12, 123–131 (2011).
  * Maxey (1987) M. R. Maxey, “The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields.” J. Fluid Mech. 174, 441–465 (1987).
  * Onishi, Baba, and Takahashi (2011) R. Onishi, Y. Baba, and K. Takahashi, “Large-scale forcing with less communication in finite-difference simulations of stationary isotropic turbulence,” J. Comput. Phys. 230, 4088–4099 (2011).
  * Obligado  _et al._ (2014) M. Obligado, T. Teitelbaum, A. Cartellier, P. Mininni, and M. Bourgoin, “Preferential concentration of heavy particles in turbulence,” J. Turbul. 15, 293–310 (2014).
  * Matsuda and Onishi (2019) K. Matsuda and R. Onishi, “Turbulent enhancement of radar reflectivity factor for polydisperse cloud droplets,” Atmos. Chem. Phys. 19, 1785 (2019).
  * Canuto  _et al._ (2007) C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang,  _Spectral methods: fundamentals in single domains_ (Springer Science & Business Media, 2007).
  * Bassenne  _et al._ (2017) M. Bassenne, J. Urzay, K. Schneider, and P. Moin, “Extraction of coherent clusters and grid adaptation in particle-laden turbulence using wavelet filters,” Phys. Rev. Fluids 2, 054301 (2017).
  * Angioni (2021) C. Angioni, “Impurity transport in tokamak plasmas, theory, modelling and comparison with experiments,” Plasma Phys. Control. Fusion 63, 073001 (2021).
  * Rumelhart, Hinton, and Williams (1986) D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by back-propagating errors,” Nature 323, 533–536 (1986).
  * Ronneberger, Fischer, and Brox (2015) O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in  _International Conference on Medical Image Computing and Computer-Assisted Intervention_ (2015) pp. 234–241.
  * Goodfellow  _et al._ (2020) I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial networks,” Commun. ACM 63, 139–144 (2020).
  * Isola  _et al._ (2017) P. Isola, J.-Y. Zhu, T. Zhou, and A. A. Efros, “Image-to-image translation with conditional adversarial networks,” in  _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ (2017) pp. 1125–1134.
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint (2014), arXiv:1412.6980 .