obtained is the $6$-point SSA with five tachyons and one tensor leg. So for
simplicity let’s first consider the five tachyons, one vector SSA
($\left(z_{1},z_{5},z_{6}\right)=\left(0,1,\infty\right)$)

$\displaystyle
A_{6}\left(\zeta_{1},k_{1};k_{2};k_{3};k_{4};k_{5};k_{6}\right)$
$\displaystyle=\int_{0}^{1}dz_{4}\int_{0}^{z_{4}}dz_{3}\int_{0}^{z_{3}}dz_{2}z_{2}^{k_{2}\cdot
k_{1}}z_{3}^{k_{3}\cdot k_{1}}z_{4}^{k_{4}\cdot
k_{1}}\left(1-z_{2}\right)^{k_{5}\cdot k_{2}}\left(1-z_{3}\right)^{k_{5}\cdot
k_{3}}$ $\displaystyle\times\left(z_{4}-z_{2}\right)^{k_{4}\cdot
k_{2}}\left(z_{4}-z_{3}\right)^{k_{4}\cdot
k_{3}}\left(z_{3}-z_{2}\right)^{k_{3}\cdot k_{2}}$
$\displaystyle\times(-\zeta_{1})\cdot\left[\frac{k_{2}}{z_{2}}+\frac{k_{3}}{z_{3}}+\frac{k_{4}}{z_{4}}+\frac{k_{5}}{1}\right]$
$\displaystyle=A_{6}^{(1)}+A_{6}^{(2)}+A_{6}^{(3)}+A_{6}^{(4)}.$ (5.69)

Mathematically, the calculation of all four terms in Eq.(5.69) are similar to
the $6$-point KN amplitude which has been proved to be a LSSA in the last
subsection. So the amplitude in Eq.(5.69) is a LSSA. According to the shifting
principle, one can generalize the above calculation and consider five
tachyons, one arbitrary tensor SSA and express them in terms of the LSSA. This
proves that the residue $R_{3,6}^{M}$ of $7$-point KN amplitude is a LSSA.
Similar consideration applies to the other residue $R_{6,3}^{M}$ of the
$7$-point KN amplitude. The other two residues $R_{4,5}^{M}$, $R_{5,4}^{M}$
contain $5$-point and $4$-point SSA which have been shown to be LSSA. We
conclude that the $7$-point KN amplitude is a LSSA. The next step is to use
the shifting principle to argue that all $7$-point SSA with arbitrary seven
tensor legs can be expressed in terms of the LSSA.

It is interesting to note that on the calculation of multi-step recursion
processes, one encounters multi-tensor lower point SSA. For example, in
addition to the $\zeta$ ($6$-point to $5$-point) in Eq.(5.23), there is
another hidden leg in $R_{3,4}^{(1),M}$ ($R_{4,3}^{(1),M}$) when one performs
the recursion from $5$-point to $4$-point. In sum, one obtains two tensor-legs
in the final $4$-point LSSA when performing $2$-step recursion
$(recursion)^{2}$ on the $6$-point KN amplitude. It is interesting to see that
the multi-tensor SSA are associated with the multi-step recursion processes.

###  V.5 N-point SSA with tensor legs

In general, we can use mathematical induction, together with the on-shell
recursion and the shifting principle, to show that all $n$-point SSA can be
expressed in terms of the LSSA. The procedure goes as following. We assume
that all $k$-point SSA ($k\leq n-1$) are LSSA, and we want to prove that all
$n$-point SSA are LSSA. To prove this, we first apply the on-shell recursion
to express the residue of $n$-point KN amplitude calculated in section III in
terms of the lower point ($k\leq n-1$) SSA which were assumed to be LSSA. See
Fig 1 above Eq.(3.27). So the $n$-point KN amplitude is a LSSA. We can then
apply the shifting principle to show that all $n$-point SSA including the
$n$-point KN amplitude are LSSA. This completes the proof.

As an example of application of the shifting principle, the $(n-1)$ tachyons,
one vector SSA can be written as
($\left(z_{1},z_{n-1},z_{n}\right)=\left(0,1,\infty\right)$)

$\displaystyle
A_{n}\left(\zeta_{1},k_{1};k_{2};k_{3};k_{4};\cdot\cdot\cdot;k_{n}\right)$
$\displaystyle=\int_{0}^{1}dz_{n-2}\int_{0}^{z_{n-2}}dz_{n-3}\cdots\int_{0}^{z_{3}}dz_{2}\text{
}\prod_{i>j=1}^{n-2}\left(z_{i}-z_{j}\right)^{k_{i}\cdot k_{j}}$
$\displaystyle\times(-\zeta_{1})\cdot\left[\frac{k_{2}}{z_{2}}+\frac{k_{3}}{z_{3}}+\frac{k_{4}}{z_{4}}+\cdot\cdot\cdot+\frac{k_{n-1}}{1}\right]$
$\displaystyle=A_{n}^{(1)}+A_{n}^{(2)}+A_{n}^{(3)}+\cdot\cdot\cdot+A_{n}^{(n-2)}$
(5.70)

where the calculation of each $A_{n}^{(j)}$ is similar to the $n$-point KN
amplitude and can be expressed in terms of LSSA.

##  VI Conclusion

In this paper, we calculate explicitly residues of all $n$-point Koba-Nielsen
(KN) amplitudes by using string on-shell recursion relation of string
scattering amplitudes (SSA). We then use mathematical induction, together with
the on-shell recursion and the shifting principle, to show that all $n$-point
SSA including the KN amplitudes can be expressed in terms of the LSSA. In
general, to explicitly calculate the LSSA form of the $n$-point KN amplitude,
one needs to apply up to ($n-4)$-step recursion $(recursion)^{n-4}$ to achieve
a $4$-point LSSA. These results demonstrates the exact $SL(K+3,\mathbb{C})$
symmetry of the tree-level open bosonic string theory.

Moreover, we derive an iteration relation among the residues of a given
$n$-point KN amplitude. This iteration relation is a generalization of the
case of simple Veneziano amplitude to the higher point SSA. We argue that this
iteration relation is related to the $SL(K+3,\mathbb{C})$ symmetry and can
presumably, in the $4$-point case, be used to soften the well-known hard
Veneziano amplitude. So we expect the $SL(K+3,\mathbb{C})$ symmetry can be
used to soften all higher $n$-point hard SSA. This mechanism is reminiscent of
symmetry principle in QFT.

In the calculation of multi-tensor SSA, we discover that the calculation of
SSA expressed in terms of the Lauricella functions is much simpler than the
traditional calculation KLT . It is also interesting to see that to express
the residues of higher point SSA in terms of the Lauricella $4$-point
function, one encounters application of the multi-step recursion processes
which are associated with the multi-tensor SSA calculation.

It will be important to extend the results of this paper to the superstring
theory. However, the known construction of the string fermion vertex involved
only the leading Regge trajectory fermion string state of the R sector PFSSA ;
Osch ; RRR . It is a nontrivial task to construct the general massive fermion
string vertex operators m1 ; m2 ; m3 ; m4 . Many questions related to the
construction of SSA involving the general massive fermion string states need
to be answered before one can further understand the symmetry of superstring
theory.

In addition to the dynamics of $4$-point SSA with higher spin particles, there
are many interesting issues related to higher point SSA. To name some, is
there any evidence of the $SL(K+3,\mathbb{C})$ symmetry other than deriving it
from SSA? In a recent publication LLYT , it was shown that the LSSA of three
tachyons and one arbitrary string states can be rederived from the deformed
cubic string field theory (SFT) Taejin . But can one use Witten SFT to obtain
the $SL(K+3,\mathbb{C})$ symmetry ? How to deal with the
$z_{4}\rightarrow\infty$ limit as was shown in Eq.(5.15) when calculating SSA
with higher spin particles putting at $z_{4}$ ? Are there linear relations
among higher point ($n\geq 5$) hard SSA similar to those of the $4$-point
case? We hope to address these issues in the next publications.

###### Acknowledgements.

We would like to thank H. Kawai and C.I. Tan for discussions which help to
clarify many issues of the LSSA. This work is supported in part by the
Ministry of Science and Technology (MoST) and S.T. Yau center of National Yang
Ming Chiao Tung University (NYCU), Taiwan.

##  Appendix A Generating function

In this appendix, we prove two relations which we used in the text to do the
residue calculation. We begin with the generating function

$e^{x\left[z+\frac{z^{2}}{2}+\frac{z^{3}}{3}+\cdots\right]}=e^{\left[xz+\frac{xz^{2}}{2}+\frac{xz^{3}}{3}+\cdots\right]}$
(A.1)

which we used in Eq.(3.1) to do the residue calculation of the Veneziano
amplitude. For our purpose here, let’s generalize this function and consider

$\displaystyle A\left(x_{1},x_{2},\cdots;z\right)$
$\displaystyle=\exp\left[B\left(x_{1},x_{2},\cdots;z\right)\right]=e^{\left[x_{1}z+\frac{x_{2}z^{2}}{2}+\frac{x_{3}z^{3}}{3}+\cdots\right]}$
$\displaystyle=\underset{m=0}{\overset{\infty}{\sum}}E_{m}\left(x_{1},x_{2},\cdots,x_{m}\right)z^{m}.$
(A.2)

Two simple applications of the functions
$E_{m}\left(x_{1},x_{2},\cdots,x_{m}\right)$ follow. If we set $x_{i}=x$ for
all $i$, we obtain

$\displaystyle
A\left(x_{1},x_{2},\cdots;z\right)=\exp\left[x\left(z+\frac{z^{2}}{2}+\frac{z^{3}}{3}+\frac{z^{4}}{4}+\frac{z^{5}}{5}\cdots\right)\right]$
$\displaystyle=E_{0}+E_{1}\left(x\right)z+E_{2}\left(x\right)z^{2}+E_{3}\left(x\right)z^{3}+E_{4}\left(x\right)z^{4}+\cdots$
(A.3)

with

$\displaystyle E_{0}$ $\displaystyle=1=\frac{\left(x\right)_{0}}{0!}$
$\displaystyle E_{1}$
$\displaystyle=\frac{1}{1}x=x=\frac{x}{1!}=\frac{\left(x\right)_{1}}{1!}$
$\displaystyle E_{2}$
$\displaystyle=\frac{1}{2}x+\frac{1}{2!}\left(\frac{1}{1}\right)^{2}x^{2}=\frac{x\left(x+1\right)}{2!}=\frac{\left(x\right)_{2}}{2!}$
$\displaystyle E_{3}$
$\displaystyle=\frac{1}{3}x+\frac{1}{1!1!}\left(\frac{1}{2}\right)\left(\frac{1}{1}\right)x^{2}+\frac{1}{3!}\left(\frac{1}{1}\right)^{3}x^{3}=\frac{x\left(x+1\right)\left(x+2\right)}{3!}=\frac{\left(x\right)_{3}}{3!}$
$\displaystyle E_{4}$
$\displaystyle=\frac{1}{4}x+\left[\frac{1}{1!1!}\left(\frac{1}{3}\right)\left(\frac{1}{1}\right)+\frac{1}{2!}\left(\frac{1}{2}\right)^{2}\right]x^{2}$
$\displaystyle+$
$\displaystyle\frac{1}{2!1!}\left(\frac{1}{1}\right)^{2}\left(\frac{1}{2}\right)x^{3}+\frac{1}{4!}\left(\frac{1}{1}\right)^{4}x^{4}=\frac{x\left(x+1\right)\left(x+2\right)\left(x+3\right)}{4!}=\frac{\left(x\right)_{4}}{4!}$
$\displaystyle\cdots$ (A.4)

The coefficient functions above find application in Eq.(3.6) which corresponds
to the residues of the Veneziano amplitude. On the other hand, if we set
$\left(x_{1},x_{2},\cdots;z\right)=\left(\alpha_{-1},\alpha_{-2},\cdots;z\right)$,
we obtain

$\displaystyle
A\left(\alpha_{-1},\alpha_{-2},\cdots;z\right)=\exp\left[\alpha_{-1}z+\frac{\alpha_{-2}z^{2}}{2}+\frac{\alpha_{-3}z^{3}}{3}+\frac{\alpha_{-4}z^{4}}{4}+\frac{\alpha_{-5}z^{5}}{5}\cdots\right]$
$\displaystyle=E_{0}+E_{1}\left(\alpha_{-1}\right)z+E_{2}\left(\alpha_{-1},\alpha_{-2}\right)z^{2}+E_{3}\left(\alpha_{-1},\alpha_{-2},\alpha_{-3}\right)z^{3}$
$\displaystyle+E_{4}\left(\alpha_{-1},\alpha_{-2},\alpha_{-3},\alpha_{-4}\right)z^{4}+...$
(A.5)

with

$\displaystyle E_{0}$ $\displaystyle=1$ $\displaystyle E_{1}$
$\displaystyle=\frac{\alpha_{-1}}{1}$ $\displaystyle E_{2}$
$\displaystyle=\frac{\alpha_{-2}}{2}+\frac{1}{2!}\left(\frac{\alpha_{-1}}{1}\right)^{2}$
$\displaystyle E_{3}$
$\displaystyle=\frac{\alpha_{-3}}{3}+\frac{1}{1!1!}\left(\frac{\alpha_{-2}}{2}\right)\left(\frac{\alpha_{-1}}{1}\right)+\frac{1}{3!}\left(\frac{\alpha_{-1}}{1}\right)^{3}$
$\displaystyle E_{4}$
$\displaystyle=\frac{\alpha_{-4}}{4}+\frac{1}{1!1!}\left(\frac{\alpha_{-3}}{3}\right)\left(\frac{\alpha_{-1}}{1}\right)+\frac{1}{2!}\left(\frac{\alpha_{-2}}{2}\right)^{2}+\frac{1}{2!1!}\left(\frac{\alpha_{-1}}{1}\right)^{2}\left(\frac{\alpha_{-2}}{2}\right)+\frac{1}{4!}\left(\frac{\alpha_{-1}}{1}\right)^{4}$
$\displaystyle\cdots$ $\displaystyle E_{M}$
$\displaystyle=\underset{\left\\{N_{t}\right\\},M=\underset{t=1}{\overset{\infty}{{\displaystyle\sum}}}tN_{t}}{\sum}\prod_{r=1}\frac{\left(\alpha_{-r}\right)^{N_{r}}}{N_{r}!r^{Nr}}$
(A.6)

where in Eq.(A.6) the sum is over partitions of $N$ into
$\left\\{N_{r}\right\\}$ with $M=\underset{r=1}{\sum}rN_{r}$ and
$N=\underset{r=1}{\sum}N_{r}$. The coefficient $E_{M}$ above reminds us of the
sum over normalized creation operators of the string spectrum at mass level
$M$.

There are two relations which we used in the text to do the residue
calculation. They are

$\displaystyle F_{M}\left(x_{1},x_{2}\right)$
$\displaystyle=\underset{b+c=M,b,c=0}{\sum}\frac{\left(T_{1}\right)_{b}}{b!}\frac{\left(T_{2}\right)_{c}}{c!}x_{1}^{b}x_{2}^{c}$
$\displaystyle=\underset{\left\\{N_{r}\right\\},M=\underset{r=1}{\sum}rN_{r}}{\sum}\prod_{r=1}\left(\frac{1}{N_{r}!r^{N_{r}}}\left[T_{1}x_{1}^{r}+T_{2}x_{2}^{r}\right]^{N_{r}}\right),$
(A.7)

and the following recurrence relation

$F_{M}\left(x_{1},x_{2}\right)=\underset{N=0}{\overset{M-1}{\sum}}F_{N}\left(x_{1},x_{2}\right)\left[T_{1}x_{1}^{M-N}+T_{2}x_{2}^{M-N}\right].$
(A.8)

To prove Eq.(A.7), we introduce the generating function

$\displaystyle G\left(x_{1},x_{2};z\right)$
$\displaystyle=\left(1-zx_{1}\right)^{-T_{1}}\left(1-zx_{2}\right)^{-T_{2}}$
$\displaystyle=\underset{b,c=0}{\sum}\frac{\left(T_{1}\right)_{b}}{b!}\frac{\left(T_{2}\right)_{c}}{c!}x_{1}^{b}x_{2}^{c}z^{b+c},$
(A.9)

and

$\displaystyle F_{M}\left(x_{1},x_{2}\right)$
$\displaystyle=\frac{1}{M!}\frac{\partial^{M}}{\partial
z^{M}}G\left(x_{1},x_{2};z\right)|_{z=0}$
$\displaystyle=\underset{b+c=M,b,c=0}{\sum}\frac{\left(T_{1}\right)_{b}}{b!}\frac{\left(T_{2}\right)_{c}}{c!}x_{1}^{b}x_{2}^{c}.$
(A.10)

On the other hand, we note that

$\displaystyle G\left(x_{1},x_{2},z\right)$
$\displaystyle=\left(1-zx_{1}\right)^{-T_{1}}\left(1-zx_{2}\right)^{-T_{2}}=\exp\left[-T_{1}\ln\left(1-zx_{1}\right)-T_{2}\ln\left(1-zx_{2}\right)\right]$
$\displaystyle=\exp\left[T_{1}\left(zx_{1}+\frac{\left(zx_{1}\right)^{2}}{2}+\frac{\left(zx_{1}\right)^{3}}{3}+\cdots\right)+T_{2}\left(zx_{2}+\frac{\left(zx_{2}\right)^{2}}{2}+\frac{\left(zx_{1}\right)^{3}}{3}+\cdots\right)\right]$
$\displaystyle=\exp\left[\left(T_{1}x_{1}+T_{2}x_{2}\right)z+\left(T_{1}x_{1}^{2}+T_{2}x_{2}^{2}\right)\frac{z^{2}}{2}+\left(T_{1}x_{1}^{3}+T_{2}x_{2}^{3}\right)\frac{z^{3}}{3}+\cdots\right]$
$\displaystyle=A\left(T_{1}x_{1}+T_{2}x_{2},T_{1}x_{1}^{2}+T_{2}x_{2}^{2},T_{1}x_{1}^{3}+T_{2}x_{2}^{3},\cdots;z\right),$
(A.11)

so

$\displaystyle F_{M}\left(x_{1},x_{2}\right)$
$\displaystyle=E_{M}\left(T_{1}x_{1}+T_{2}x_{2},T_{1}x_{1}^{2}+T_{2}x_{2}^{2},\cdots,T_{1}x_{1}^{M}+T_{2}x_{2}^{M}\right)$
$\displaystyle=\sum_{\left\\{N_{r}\right\\}\text{,with
}M=\underset{r}{\sum}rN_{r}}\prod_{r=1}\left(\frac{\left(T_{1}x_{1}^{r}+T_{2}x_{2}^{r}\right)^{N_{r}}}{N_{r}!r^{N_{r}}}\right)$
(A.12)

where we have used Eq.(A.6) to obtain Eq.(A.12). This completes the proof of
Eq.(A.7). It is straightforward to generalize Eq.(A.7) to

$\displaystyle F_{M}\left(x_{1},\cdots,x_{n}\right)$
$\displaystyle=\underset{\begin{array}[c]{c}a_{1}+\cdots+a_{n}=M\\\
a_{1},\cdots,a_{n}=0\end{array}}{\sum}\frac{\left(T_{1}\right)_{a_{1}}}{a_{1}!}\cdots\frac{\left(T_{n}\right)_{a_{n}}}{a_{n}!}x_{1}^{a_{1}}\cdots
x_{n}^{a_{n}}$ (A.15)
$\displaystyle=\underset{\left\\{N_{r}\right\\}\text{,with fixed
}M=\underset{r}{\sum}rN_{r}}{\sum}\prod_{r=1}\left(\frac{1}{N_{r}!r^{N_{r}}}\left[T_{1}x_{1}^{r}+\cdots+T_{n}x_{n}^{r}\right]^{N_{r}}\right).$
(A.16)

It can be shown that $E_{m}\left(x_{1},\cdots,x_{m}\right)$ satisfies a
iteration relation

$\displaystyle E_{m}\left(x_{1},\cdots,x_{m}\right)$
$\displaystyle=\frac{1}{m!}\frac{\partial^{m}}{\partial
z^{m}}A\left(x_{1},\cdots;z\right)|_{z=0}=\frac{1}{m!}\frac{\partial^{m-1}}{\partial
z^{m-1}}\left(\frac{\partial B}{\partial z}e^{B}\right)|_{z=0}$
$\displaystyle=\frac{1}{m}\underset{n=0}{\overset{m-1}{\sum}}\frac{1}{n!}\left(\frac{\partial^{n}}{\partial
z^{n}}e^{B}\right)\frac{1}{\left(m-n-1\right)!}\frac{\partial^{m-1-n}}{\partial
z^{m-1-n}}\left(\frac{\partial B}{\partial z}\right)|_{z=0}$
$\displaystyle=\frac{1}{m}\underset{n=0}{\overset{m-1}{\sum}}\frac{1}{n!}\left(\frac{\partial^{n}}{\partial
z^{n}}e^{B}\right)|_{z=0}\left[x_{m-n}+O\left(z\right)\right]|_{z=0}=\frac{1}{m}\underset{n=0}{\overset{m-1}{\sum}}E_{n}\left(x_{1},\cdots,x_{n}\right)x_{m-n}.$
(A.17)

If we put $x_{i}=\alpha_{-i}$ for all $i$, we obtain

$E_{m}\left(\alpha_{-1},\alpha_{-2},\cdots;\alpha_{-m}\right)=\frac{1}{m}\underset{n=1}{\overset{m-1}{\sum}}E_{n}\left(\alpha_{-1},\alpha_{-2},\cdots;\alpha_{-n}\right)\alpha_{-\left(m-n\right)}.$
(A.18)

Similarly, $F_{M}\left(x_{1},x_{2}\right)$ obeys a iteration relation

$\displaystyle
F_{M}\left(x_{1},x_{2}\right)=\frac{1}{M!}\frac{\partial^{M}}{\partial
z^{M}}G\left(x_{1},x_{2};z\right)|_{z=0}$
$\displaystyle=\frac{1}{M!}\frac{\partial^{M}}{\partial
z^{M}}\left(1-zx_{1}\right)^{-T_{1}}\left(1-zx_{2}\right)^{-T_{2}}|_{z=0}$
$\displaystyle=\frac{1}{M!}\frac{\partial^{M-1}}{\partial
z^{M-1}}\left[\left(1-zx_{1}\right)^{-T_{1}}\left(1-zx_{2}\right)^{-T_{2}}\left(\frac{x_{1}T_{1}}{1-zx_{1}}+\frac{x_{1}T_{2}}{1-zx_{2}}\right)\right]|_{z=0},$
(A.19)

which can be written as

$\displaystyle=\frac{1}{M!}\underset{N=0}{\overset{M-1}{\sum}}\left[\begin{array}[c]{c}\frac{\left(M-1\right)!}{N!\left(M-1-N\right)!}\frac{\partial^{N}}{\partial
z^{N}}\left(\left(1-zx_{1}\right)^{-T_{1}}\left(1-zx_{2}\right)^{-T_{2}}\right)|_{z=0}\\\
\times\frac{\partial^{M-1-N}}{\partial
z^{M-1-N}}\left(\frac{x_{1}T_{1}}{1-zx_{1}}+\frac{x_{1}T_{2}}{1-zx_{2}}\right)|_{z=0}\end{array}\right]$
(A.22)
$\displaystyle=\frac{1}{M}\underset{N=0}{\overset{M-1}{\sum}}\frac{1}{N!}\frac{\partial^{N}}{\partial
z^{N}}\left(\left(1-zx_{1}\right)^{-T_{1}}\left(1-zx_{2}\right)^{-T_{2}}\right)|_{z=0}\left(T_{1}x_{1}^{M-N}+T_{2}x_{2}^{M-N}\right)$
$\displaystyle=\frac{1}{M}\underset{N=0}{\overset{M-1}{\sum}}F_{N}\left(x_{1},x_{2}\right)\left(T_{1}x_{1}^{M-N}+T_{2}x_{2}^{M-N}\right).$
(A.23)

This completes the proof of Eq.(A.8). It is straightforward to generalize
Eq.(A.8) to

$F_{M}\left(x_{1},x_{2},\cdots,x_{n}\right)=\underset{N=0}{\overset{M-1}{\sum}}F_{N}\left(x_{1},x_{2},\cdots,x_{n}\right)\left[T_{1}x_{1}^{M-N}+T_{2}x_{2}^{M-N}+\cdots+T_{n}x_{n}^{M-N}\right].$
(A.24)

## References

  * (1) MB Green, JH Schwarz and E Witten. Superstring theory, v. 1. Cambridge University,Cambridge, (1987).
  * (2) Gregory Moore.  Finite in all directions.  arXiv preprint hep-th/9305139, 1993.
  * (3) Gregory Moore.  Symmetries of the bosonic string S-matrix.  arXiv preprint hep-th/9310026, 1993.
  * (4) Chuan-Tsung Chan, Shoichi Kawamoto, and Dan Tomino.  To see symmetry in a forest of trees.  Nucl. Phys. B, 885:225–265, 2014.
  * (5) Sheng-Hong Lai, Jen-Chi Lee, and Yi Yang.  The lauricella functions and exact string scattering amplitudes.  Journal of High Energy Physics, 2016(11):62, 2016.
  * (6) S. H. Lai, J. C. Lee and Yi Yang, ”The $SL(K+3,C)$ Symmetry of the Bosonic String Scattering Amplitudes”, Nucl. Phys. B 941 (2019) 53-71.
  * (7) Willard Miller Jr.  Symmetry and separation of variables. Addison-Wesley, Reading, Massachusetts,  1977\.
  * (8) S. H. Lai, J. C. Lee and Yi Yang, ”Solving Lauricella String Scattering Amplitudes through recurrence relations”, JHEP 09 (2017) 130.
  * (9) S. H. Lai, J. C. Lee and Yi Yang, ”Recent developments of the Lauricella string scattering amplitudes and their exact $SL(K+3,C)$ Symmetry”, Symmetry 13 (2021) 454, arXiv:2012.14726 [hep-th].
  * (10) Steven B Giddings, David J Gross, and Anshuman Maharana.  Gravitational effects in ultrahigh-energy string scattering.  Physical Review D, 77(4):046001, 2008.
  * (11) Daniele Amati, M Ciafaloni, and G Veneziano.  Superstring collisions at Planckian energies.  Phys. Lett. B, 197(1):81–88, 1987.
  * (12) Daniele Amati, M Ciafaloni, and G Veneziano.  Classical and quantum gravity effects from Planckian energy superstring collisions.  Int. J. of Mod. Phys. A, 3(07):1615–1661, 1988.
  * (13) Daniele Amati, Marcello Ciafaloni, and Gabriele Veneziano.  Can spacetime be probed below the string size?  Phys. Lett. B, 216(1):41–47, 1989.
  * (14) M Soldate.  Partial-wave unitarity and closed-string amplitudes.  Phys. Lett. B, 186(3):321–327, 1987.
  * (15) IJ Muzinich and M Soldate.  High-energy unitarity of gravitation and strings.  Phys. Rev. D, 37(2):359, 1988.
  * (16) Richard C Brower, Joseph Polchinski, Matthew J Strassler, and Chung-I Tan.  The Pomeron and gauge/string duality.  JHEP, 2007(12):005, 2007.
  * (17) Richard C Brower, Horatiu Nastase, Howard J Schnitzer, and Chung-I Tan.  Analyticity for multi-Regge limits of the Bern–Dixon–Smirnov amplitudes.  Nucl. Phys. B, 822(1):301–347, 2009.
  * (18) David J Gross and Paul F Mende.  The high-energy behavior of string scattering amplitudes.  Phys. Lett. B, 197(1):129–134, 1987.
  * (19) David J Gross and Paul F Mende.  String theory beyond the Planck scale.  Nucl. Phys. B, 303(3):407–454, 1988.
  * (20) David J Gross and JL Manes.  The high energy behavior of open string scattering.  Nucl. Phys. B, 326(1):73–107, 1989.
  * (21) David J. Gross.  High-Energy Symmetries of String Theory.  Phys. Rev. Lett., 60:1229, 1988.
  * (22) D. J. Gross.  Strings at superPlanckian energies: In search of the string symmetry.  In _In *London 1988, Proceedings, Physics and mathematics of strings* 83-95. (Philos. Trans. R. Soc. London A329 (1989) 401-413)._ , 1988.
  * (23) Chuan-Tsung Chan and Jen-Chi Lee.  Stringy symmetries and their high-energy limits.  Phys. Lett. B, 611(1):193–198, 2005.
  * (24) Chuan-Tsung Chan and Jen-Chi Lee.  Zero-norm states and high-energy symmetries of string theory.  Nucl. Phys. B, 690(1):3–20, 2004.
  * (25) Chuan-Tsung Chan, Pei-Ming Ho, and Jen-Chi Lee.  Ward identities and high energy scattering amplitudes in string theory.  Nucl. Phys. B, 708(1):99–114, 2005.
  * (26) Chuan-Tsung Chan, Pei-Ming Ho, Jen-Chi Lee, Shunsuke Teraguchi, and Yi Yang.  High-energy zero-norm states and symmetries of string theory.  Phys. Rev. Lett., 96(17):171601, 2006.
  * (27) Chuan-Tsung Chan, Pei-Ming Ho, Jen-Chi Lee, Shunsuke Teraguchi, and Yi Yang.  Solving all 4-point correlation functions for bosonic open string theory in the high-energy limit.  Nucl. Phys. B, 725(1):352–382, 2005.
  * (28) Jen-Chi Lee.  New symmetries of higher spin states in string theory.  Physics Letters B, 241(3):336–342, 1990.
  * (29) Jen-Chi Lee and Burt A Ovrut.  Zero-norm states and enlarged gauge symmetries of the closed bosonic string with massive background fields.  Nucl. Phys. B, 336(2):222–244, 1990.
  * (30) Jen-Chi Lee.  Decoupling of degenerate positive-norm states in string theory.  Phys. Rev. Lett., 64(14):1636, 1990.
  * (31) Jen-Chi Lee and Yi Yang.  Review on high energy string scattering amplitudes and symmetries of string theory.  arXiv preprint arXiv:1510.03297, 2015.
  * (32) Jen-Chi Lee and Yi Yang.  Overview of high energy string scattering amplitudes and symmetries of string theory.  Symmetry, 11(8):1045, 2019.
  * (33) Bagchi, Arjun, Aritra Banerjee, and Shankhadeep Chakrabortty. ”Rindler Physics on the String Worldsheet.” Physical Review Letters 126.3 (2021): 031601.
  * (34) Eberhardt, Lorenz. ”Partition functions of the tensionless string.” arXiv preprint arXiv:2008.07533 (2020).
  * (35) Bagchi, Arjun, et al. ”A tale of three—tensionless strings and vacuum structure.” Journal of High Energy Physics 2020.4 (2020): 1-53.
  * (36) Lee, Seung-Joo, Wolfgang Lerche, and Timo Weigand. ”Tensionless strings and the weak gravity conjecture.” Journal of High Energy Physics 2018.10 (2018): 1-83.
  * (37) Bagchi, Arjun, et al. ”Inhomogeneous tensionless superstrings.” Journal of High Energy Physics 2018.2 (2018): 1-33.
  * (38) Yu, Ming, Chi Zhang, and Yao-Zhong Zhang. ”One loop amplitude from null string.” Journal of High Energy Physics 2017.6 (2017): 1-18.
  * (39) Hohm, Olaf, Usman Naseer, and Barton Zwiebach. ”On the curious spectrum of duality invariant higher-derivative gravity.” Journal of High Energy Physics 2016.8 (2016): 1-31.
  * (40) Bagchi, Arjun, Shankhadeep Chakrabortty, and Pulastya Parekh. ”Tensionless superstrings: view from the worldsheet.” Journal of High Energy Physics 2016.10 (2016): 1-24.
  * (41) Bagchi, Arjun, Shankhadeep Chakrabortty, and Pulastya Parekh. ”Tensionless strings from worldsheet symmetries.” Journal of High Energy Physics 2016.1 (2016): 158.
  * (42) Gaberdiel, Matthias R., and Rajesh Gopakumar. ”Higher spins & strings.” Journal of High Energy Physics 2014.11 (2014): 44.
  * (43) Sagnotti, A., and M. Tsulaia. ”On higher spins and the tensionless limit of string theory.” Nuclear Physics B 682.1-2 (2004): 83-116.
  * (44) Yung-Yeh Chang, Bo Feng, Chih-Hao Fu, Jen-Chi Lee,Yihong Wang and Yi Yang, ”A note on on-shell recursion relation of string amplitudes”, JHEP 02 (2013) 028.
  * (45) R. H. Boels, D. Marmiroli and N. A. Obers, “On-shell Recursion in String Theory,” JHEP 1010 (2010) 034 [arXiv:1002.5029 [hep-th]].
  * (46) C. Cheung, D. O’Connell and B. Wecht, “BCFW Recursion Relations and String Theory,” JHEP 1009 (2010) 052 [arXiv:1002.4674 [hep-th]].
  * (47) R. Britto, F. Cachazo and B. Feng, “New recursion relations for tree amplitudes of gluons,” Nucl. Phys. B 715 (2005) 499 [hep-th/0412308].
  * (48) R. Britto, F. Cachazo, B. Feng and E. Witten, “Direct proof of tree-level recursion relation in Yang-Mills theory,” Phys. Rev. Lett. 94 (2005) 181602 [hep-th/0501052].
  * (49) Sheng-Hong Lai, Jen-Chi Lee, Taejin Lee and Yi Yang, Phys. Lett. B 776 (2018) 150-157.
  * (50) Taejin Lee, Phys. Lett. B 768 (2017) 248.
  * (51) J AppellPandKampé de Fériet.  Fonctions hypergéométriques et hypersphériques, 1926.
  * (52) Ko, S.-L.; Lee, J.-C.; Yang, Y.  Patterns of high energy massive string scatterings in the regge regime.  J. High Energy Phys. 2009, _2009_ , 28.
  * (53) Lee, J.-C.; Yan, C.H.; Yang, Y.  High-energy string scattering amplitudes and signless Stirling number identity.  SIGMA Symmetry Integr. Geom. Methods Appl. 2012, _8_ , 45.
  * (54) H. Kawai, David C Lewellen, and S-HH Tye.  A relation between tree amplitudes of closed and open strings.  Nucl. Phys. B, 269(1):1–23, 1986.
  * (55) Sheng-Hong Lai, Jen-Chi Lee, and Yi Yang, work in progress.
  * (56) S. H. Lai, J. C. Lee and Yi Yang, Phys. Lett. B797, 134812 (2019).
  * (57) O. Schlotterer, ”Scattering amplitudes in open superstring theory”, Fortschr. Phys.60,No.5, 373-691 (2012).
  * (58) I.G. Koh, W. Troost and A. Van Proeyen, Nucl. Phys. B292, 201 (1987).
  * (59) M. Bianchi, L. Lopez, and R. Richter, ”On stable higher spin states in Heterotic String Theories,” JHEP 1103 (2011) 051
  * (60) A. Hanany, D. Forcella, and J. Troost, ”The Covariant perturbative string spectrum”, Nucl.Phys. B846 (2011) 212.
  * (61) Wan-Zhe Feng, Dieter Lust, Oliver Schlotterer, Stephan Stieberger, Tomasz R. Taylor, Nucl.Phys.B843, 570 (2011).
  * (62) M. Vasiliev, ”Nonlinear equations for symmetric massless higher spin fields in (A)dS(d)”, Phys.Lett. B567 (2003) 139.


# Data-driven model for divertor plasma detachment prediction

Ben Zhu<EMAIL_ADDRESS>Menglong Zhao  Harsh Bhatia  Xue-qiao Xu  Peer-Timo
Bremer  William Meyer  Nami Li  Thomas Rognlien  Lawrence Livermore National
Laboratory, Livermore, CA 94550, USA

###### Abstract

We present a fast and accurate data-driven surrogate model for divertor plasma
detachment prediction leveraging the latent feature space concept in machine
learning research. Our approach involves constructing and training two neural
networks. An autoencoder that finds a proper latent space representation (LSR)
of plasma state by compressing the multi-modal diagnostic measurements, and a
forward model using multi-layer perception (MLP) that projects a set of plasma
control parameters to its corresponding LSR. By combining the forward model
and the decoder network from autoencoder, this new data-driven surrogate model
is able to predict a consistent set of diagnostic measurements based on a few
plasma control parameters. In order to ensure that the crucial detachment
physics is correctly captured, highly efficient 1D UEDGE model is used to
generate training and validation data in this study. Benchmark between the
data-driven surrogate model and UEDGE simulations shows that our surrogate
model is capable to provide accurate detachment prediction (usually within a
few percent relative error margin) but with at least four orders of magnitude
speed-up, indicating that performance-wise, it has the potential to facilitate
integrated tokamak design and plasma control. Comparing to the widely used
two-point model and/or two-point model formatting, the new data-driven model
features additional detachment front prediction and can be easily extended to
incorporate richer physics. This study demonstrates that the complicated
divertor and scrape-off-layer plasma state has a low-dimensional
representation in latent space. Understanding plasma dynamics in latent space
and utilizing this knowledge could open a new path for plasma control in
magnetic fusion energy research.

##  1 Introduction

One of the critical concerns for a magnetic fusion power plant is the intense
plasma heat flux $q_{\parallel}$ that flows very rapidly along magnetic field
lines to localized material divertor plates. For unmitigated conditions, this
heat flux can exceed $10~{}GW/m^{2}$(Kuang et al., 2020) based on
extrapolations from current experimental databases and challenges even the
most highly developed modern materials with thermo-mechanical engineering
limits at around $10~{}MW/m^{2}$. So far, the most successful method of
reducing this peak heat flux in tokamaks is achieved by detachment \- that is
to use neutral gas to buffer the incoming plasma heat flux and radiate much of
the incoming power over a much wider area of the surrounding walls as the
recombined neutrals and the radiation photons freely stream across the
magnetic field rather than being focused to the plate as for the plasma. In
other words, upstream high temperature plasma is extinguished by high density
neutral gas before it strikes the divertor plates; hence, it is detached.

Because detachment phenomenon involves nonlinear coupled atomic, molecular and
plasma physics, it is challenging to obtain accurate and fast detachment
predictions (i.e., predict divertor plasma quantities such as density and
temperature for prescribed upstream plasma condition). On one hand, theory-
based models with simplifications such as the semi-analytical two-point model
formatting (2PMF) (Stangeby, 2018) derives analytical formula to calculate
divertor plasma density and temperature. These over-simplified 0D models are
fast and valuable in terms of providing physics insights, but the
corresponding prediction is somewhat crude and relies on ad-hoc or data-fitted
momentum and energy loss fraction coefficients. On the other hand, sophistic
tokamak edge transport codes, such as 2D UEDGE (Rognlien et al., 1994) and
SOLPS-ITER (Wiesen et al., 2015) that incorporate a higher level of complexity
and physical effects could give fairly accurate predictions but require a few
days to a few months to finish one run. Of course, with a large enough neutral
gas and/or impurity injection, detachment almost always occurs; however, “too
much” detachment can lead to disruption – a rapid termination of plasma
confinement, which could also damage the machine. Thus, predicting,
establishing and maintaining a proper degree of plasma detachment, along with
a stable and properly positioned “detachment front”, is a non-trivial
requirement. To make this challenge ever greater, current detachment control
in existing tokamaks is accomplished only with the aid of in situ diagnostics
such as infrared camera (Maurizio et al., 2017) or Langmuir probe (Eldon et
al., 2022) that would not survive in a reactor environment despite that the
success of future magnetic fusion reactors, such as ITER, rely on detachment
to mitigate excessive heat load on divertor targets. Therefore, attaining a
reliable and fast detachment prediction model and control algorithm is
imperative to the magnetic fusion energy research.

To address this challenge, we propose to build a detachment prediction model
with the data-driven machine learning approach. For the last few years,
machine learning methods have gained a lot of interest in magnetized plasma
research, as summarized in latest review paper by Anirudh et al. (2022).
Various neural networks have been explored and applied to assorted
applications of experimental and computational fusion plasma study, e.g., to
accelerate data processing (van den Berg et al., 2018), to interpret critical
plasma parameter directly from diagnostics (Samuell et al., 2021), to aid
plasma shaping via magnetic field control (Degrave et al., 2022), to forecast
disruption (Kates-Harbeck et al., 2019; Montes et al., 2019; Rea et al.,
2020), to speed-up nonlinear Fokker-Planck-Landau collision operator in
particle simulations (Miller et al., 2021), to apply closures in fluid models
(Ma et al., 2020; Maulik et al., 2020; Wang et al., 2020), to discover
governing physics-based partial observations (Mathews et al., 2021), and so
on. The growing applications of machine learning techniques on plasma physics
indicates that data-driven science could provide an alternative approach to
scientific challenges.

In this paper, we report a fast yet accurate model-based data-driven
detachment prediction model. The rest of the paper is organized as follows.
The methodology of our approach is discussed in Section 2. Section 3 describes
the 1D UEDGE model used in the study and outlines the data generation process.
Section 4 presents the architecture of neural networks and hyper-parameters
used in this study. The results of our trained models are disseminated in
Section 5. Section 6 compares 1D UEDGE simulations and our predictive model
with the basic two-point model and two-point model formatting. Application of
our model, its limitations and future work are briefly discussed in section 7.
Finally, section 8 summarizes the paper.

##  2 Methodology

Figure 1: Three descriptions of plasma state.

Steady-state plasma states are described in two different ways: (1) by control
parameters ${\mathbf{x}}$, or (2) by diagnostic measurements ${\mathbf{y}}$.
For examples, if we neglect any bifurcation or hysteresis phenomena for now, a
tokamak plasma is determined by a complete set of engineering and discharge
parameters or can be characterized with of diagnosed plasma quantities in
detail. Analogously, in numerical simulations, a set of model inputs with
proper boundary conditions yields a plasma state that can also be diagnosed
synthetically. A large portion of magnetized plasma research in fact is to
discover and exploit the relation between these two descriptions, i.e.,
predicting laboratory and fusion plasma behavior under certain discharge
parameters or finding an optimal control setting for a desired plasma state.
In most situations, this is not a trivial task because of the strongly
nonlinear nature of magnetized plasmas. On one hand, sophisticated numerical
models are developed to better understand the physics which governs the
complex plasma system, and hopefully to achieve true first principle-based
predictive capability. On the other hand, empirical and simplified models are
proposed to provide tools for daily machine operations, experiment planning
and future device design. Often there is a gap in between them - empirical and
simplified models are fast but not necessarily accurate, while sophisticated
numerical models can provide more reliable predictions but too slow to be used
for practical application purpose. A fast yet accurate surrogate model can
bridge this gap. In this paper, we develop such a surrogate model for divertor
plasma detachment prediction with data-driven approach.

Unlike most data-driven surrogate models that directly connect two states
${\mathbf{x}}$ and ${\mathbf{y}}$, we take the indirect approach proposed by
(Anirudh et al., 2020). As shown in Figure 1, in addition to two conventional
descriptions of plasma, here we propose a third description — that in a latent
space — termed as latent space representation (LSR). Latent space, sometimes
also referred to as feature space, is a widely used concept in machine
learning research in which items with similar features are positioned closely.
We first find the LSR, ${\mathbf{z}}$, of plasma by compressing diagnostics
(e.g., synthetic diagnostics mimicking physical diagnostics, such as Langmuir
probe, Thomson scattering, bolometer/radiation measurement, etc.) through an
autoencoder (described in 4.1); then forward and inverse models can be trained
to relate control parameters ${\mathbf{x}}$ (engineering or numerical model
input parameters such as heating power, gas puffing rate/upstream density,
etc) and ${\mathbf{z}}$. Hence, instead of constructing neural networks
${\mathbf{x}}\to{\mathbf{z}}$ or ${\mathbf{z}}\to{\mathbf{x}}$, our forward
and inverse predictions are ${\mathbf{x}}\to{\mathbf{z}}\to{\mathbf{y}}$ and
${\mathbf{y}}\to{\mathbf{z}}\to{\mathbf{x}}$. Comparing to the direct
approach, the study (Anirudh et al., 2020) shows that indirect approach has
several advantages, including improved predictive performance, and more data
efficient. In this paper, we focus on forward prediction (i.e.,
${\mathbf{x}}\to{\mathbf{z}}\to{\mathbf{y}}$), whereas the inverse prediction
(optimization problem) is subject to future work.

##  3 Physics model and data

In this section, we briefly describe the tokamak edge model used to simulate
divertor detachment and outline the data generation process.

###  3.1 1D UEDGE model

UEDGE (Rognlien et al. (1994); Rognlien & Rensink (2002)) is a finite-volume
simulation code to capture plasma and neutral transport. UEDGE has been used
extensively to model the edge region of present-day and future tokamak devices
for almost three decades. It has been utilized to study the divertor
detachment physics on many machines, such as Alcator C-Mod (Wising et al.
(1997)), DIII-D (Porter et al. (1996)), ITER (Wising et al. (1996)), SPARC
(Ballinger et al. (2021)). UEDGE solves multi-species fluid transport model
for plasma and neutral gas, including important atomic physics, such as line
radiation from ionization, excitation, and recombination processes in
realistic tokamak equilibria. With the fully implicit Newton-Krylov scheme and
adaptive stepping to advance all model equations in time, UEDGE is able to
quickly evolve the system to a steady state.

In this proof-of-principle study, highly flexible and efficient 1D UEDGE code
is used to quickly generate a large training data set that covers a wide range
of parameter space beyond normal tokamak operation. Unlike the commonly used
2D UEDGE model, which simulates the entire tokamak periphery, including closed
flux surface region, scrape-off-layer (SOL), and private flux surface region,
the 1D UEDGE in this study models only one flux-tube in the low-field side
SOL, starting from the outer mid-plane to the outer divertor target plate, as
depicted in Figure 2(a). Nevertheless, 1D UEDGE code solves the same set of
equations as its 2D counterpart except radial transport and drifts. Plasmas
are allowed to transport in the parallel direction that is along the magnetic
field whereas neutrals transport poloidally along the flux-tube. As a result,
the essence of plasma detachment physics is retained despite the one flux-tube
assumption. In fact, simplified 1D geometry disentangles many secondary
effects and is sometimes preferred to better quantify the role of different
physics processes played in the detachment (Dudson et al. (2019)). In 1D UEDGE
setup, spatial variation of the magnetic field $\boldsymbol{B}$ is retained,
i.e. the flux expansion along a flux tube in the SOL is considered in the 1D
UEDGE mesh. Power, or the energy influx, across the core boundary is injected
at the outer midplane (upstream), and is equally partitioned between electrons
and ions. Plasma particles are also assumed to be fueled at upstream. As
plasma moves downside the flux tube, it recombines to neutrals which
eventually are recycled at the divertor target plate. In this study, the
particle and momentum recycling coefficients are set to be $0.99$.

Figure 3 shows the typical attached (in blue) and detached (in orange) plasma
in 1D UEDGE simulations. Compared with attached plasma, detached plasma
normally has a lower electron temperature $T_{e}$, a higher neutral gas
density $n_{g}$ at the divertor target, and a distinct detachment front, or
radiation peak away from the divertor target. Not surprisingly, 1D UEDGE code
is able to qualitatively reproduce the target ion saturation current density
$J_{\parallel}^{\text{sat}}$ rollover (e.g., Figure 2 in Loarte et al. (1998))
and electron temperature cliff (e.g., Figure 5(a) in McLean et al. (2015))
phenomena observed in detachment experiments as shown in Figure 4. We admit
that more sophisticated models such as 2D UEDGE and SOLPS-ITER considering
additional details in principle yield more trustworthy results; but the
fundamental physics picture of plasma detachment - that is plasma radiates
enough power to reach a low-enough temperature (a few $eV$) and to recombine
into neutral gas will not change. Rather than developing a more comprehensive
data-driven model with training data from state-of-the-art numerical
simulations, here we focus on exploring the data-driven approach for
detachment prediction with a somewhat simplified physics-based model (1D
UEDGE) that could cover lots of plasma discharge scenarios relatively cheap
computationally.

Figure 2: (a) 2D (black) v.s. 1D (red) UEDGE simulation meshes, and (b) sketch
of 1D UEDGE simulation setup. Figure 3: Examples of (a) plasma density, (b)
electron temperature, (c) neutral density, and (d) radiation profile of
attached (blue) and detached (orange) divertor plasma from 1D UEDGE
simulations.  Figure 4: Reproductions of (a) $J_{\parallel}^{\text{sat}}$
rollover, and (b) $T_{e}$ cliff features in 1D UEDGE density scan simulations.

###  3.2 Data generation

Although there are more than one dozen independent geometry and plasma
parameters can be adjusted in 1D UEDGE model, here we only vary three most
relevant parameters, namely, upstream density $n_{e,u}$, injection power
$P_{\text{inj}}$ and impurity fraction $f_{Z}$. These three parameters
$(n_{e,u},P_{\text{inj}},f_{z}$) will be referred as “model inputs”
thereafter.

In this study, we use DIII-D tokamak geometry with a fixed divertor leg length
$L_{\text{leg}}$ $=0.2696~{}m$, defined as the poloidal distance between the
$X$-point to the outer divertor target as shown in Figure2(b). We uniformly
sample $n_{e,u}\in[1,7]\times 10^{19}~{}m^{-3}$, $P_{\text{inj}}$
$\in[1,10]~{}MW$ and $f_{Z}\in 0-10\%$. Because the first wall material of
DIII-D tokamak is graphite, here the impurity species is set to be carbon.
With 60 sample points for $n_{e,u}$ and $f_{Z}$, and 40 sample points for
$P_{\text{inj}}$, 144,000 cases in total are launched. However, not all of
these cases can reach a physically feasible steady state solution within a
limited simulation time. Hence, only 111,598 converged simulation cases with
good power balance relation are accepted as the training data set. Similarly,
an independent validation data set with a coarser sampling rate ($50\times
50\times 30$ and 57,655 accepted solutions) is generated. Unless stated
otherwise, the neural networks presented in this paper are trained and
validated based on the training data set, only the combined forward prediction
model are validated with the validation data set.

data set | $n_{\text{total}}$ | $n_{\text{detach}}(\%)$ | $n_{\text{attach}}(\%)$
---|---|---|---
training | 111,598 | 56,047(50.22) | 55,551(49.78)
validation | 57,655 | 28,468(49.38) | 29,187(50.62)
Table 1: Numbers of detached and attached cases and their percentages in the
training and validation data sets.

To avoid manual intervention of hundreds of thousands simulations
individually, all the UEDGE runs are launched and managed by Merlin (Peterson
et al., 2022) \- a machine learning workflow framework designed for large
scale HPC environment. It takes 4 days with 200 CPUs to get about $170,000$
converged cases. The data generation time could be further shorten if more
CPUs are utilized.

For the diagnostic description of plasma state in this 1D setup, five
synthetic measurements are taken across the simulation domain, including the
upstream electron temperature $T_{e,u}$, ion saturation current density
$J_{\parallel}^{\text{sat}}$, electron density $n_{e,t}$ and temperature
$T_{e,t}$ at divertor target, as well as the radiation profile $P_{rad}$. As
illustrated in Figure 2 (b), $T_{e,u}$ mimics the mid-plane Thompson
scattering measurement, $J_{\parallel}^{\text{sat}}$ mimics the Langmuir probe
measurement, $n_{e,t}$ and $T_{e,t}$ come from either divertor Thompson
scattering or Langmuir probe, and $P_{rad}$ is deduced from either bolometer
or C-III emission. Note here $J_{\parallel}^{\text{sat}}$, $T_{e,u}$,
$n_{e,t}$ and $T_{e,t}$ are scalars, while $P_{rad}$ is a 1D profile. These
diagnostic measurements are collected for all 111,598 cases and then
normalized for model training.

Figure 5: Histogram and cumulative distribution of electron temperature at
divertor target $T_{e,t}$ from the training data set (111,598 cases in total).

Prior model training, all cases are labeled as either “detached” or “attached”
because there is no “partial detached” state in the 1D system. Here the choice
of detachment criterion is straightforward. As shown in Figure 5, electron
temperature at the divertor target $T_{e,t}$ has two distinct distributions
around $T_{e,t}=2.1~{}eV$. The gap at $2.1~{}eV$ is caused by the temperature
cliff phenomenon showed in Figure 4(b). Therefore, 56,047 cases with
$T_{e,t}<2.1~{}eV$ are labeled “detached”; while the rest 55,551 cases are
“attached” cases. As shown in Table 1, nearly half of the cases in the
validation data set are labeled “detached”, similar to the percentage of
“detached” cases in the training set. In DIII-D experiment, the temperature
cliff occurs around $5~{}eV$, this apparent discrepancy is mainly due to our
simplified 1D setup. The temperature cliff in our simulations is consistent
with the sonic transition point where the ion parallel velocity reaches sonic
moving away from the target. A thorough investigation of temperature cliff in
1D is beyond the scope of this paper and will be presented in a separated
paper by Zhao et al. (2022). We would like to point out that labeling cases is
primarily for better quantify models’ accuracy in Section 5; it doesn’t affect
the training process at all. Also, the lower $T_{e,t}\sim 0.22~{}eV$ limit
shown in Figure 5 is likely related to the atomic physics. When
$T_{e}<1~{}eV$, the ionization rate decreases while the recombination rate
increases dramatically. In this temperature range, electron collisionality is
very high such that three-body recombination becomes dominant. This process
moves plasma potential energy to electron thermal energy, stopping $T_{e}$
decreasing further.

##  4 Development of Machine Learning Models

The goal of our ML-based approach is to predict the outputs of UEDGE
simulations directly from a set of given model inputs, without actually
performing UEDGE simulations and with significant speed up. As has been shown
in literature Anirudh et al. (2020), this task becomes substantially easier
and more tractable if the prediction is made through an intermediate and
reduced dimensional space also learned through ML (see Figure 1 and Figure 6).
Therefore, in this work, we utilize two ML models. The first ML model — a
_variational autoencoder_ — is used to generate an invertible mapping from ODs
to a reduced-dimensional _“latent space”_ , whereas the second ML model — a
_surrogate model_ — is used to predict UEDGE outputs in this latent space,
given the inputs to UEDGE.

Figure 6: Sketches of (a) $\beta$-variational autoencoder, (b) forward
multilayer perceptron (MLP), and (c) combined data-driven model.

###  4.1 Latent Space Identification using Autoencoders

The basic idea behind latent space identification (also referred to as
dimensionality reduction or data compression) is straightforward. Given an
input data point, the goal of an autoencoder (AE) is to reduce it to fewer
dimensions such that the reduced representation allows reconstructing the
original input as faithfully as possible. An AE comprises two parts: (1) an
encoder that learns to meaningfully encode the data into a set of values — a
compact representation of the input (i.e., latent coordinates), and (2) a
decoder that attempts to reconstruct the original inputs from this latent
space encoding. By forcing the reconstruction of data, an AE attempts to learn
a suitable latent space that captures all necessary degrees of freedom while
discarding trivial variations, noise, and redundant correlations in the data.

In this work, we utilize two types of autoencoders to create a suitable latent
space representation in two steps. First, we develop a (vanilla) AE to
identify an appropriate level of reduction, i.e., the dimensionality of the
latent space, and a suitable neural network architecture. Next, we turn this
autoencoder into a special type of autoencoder, called the $\beta$-variational
autoencoder (VAE), to create the final latent space representation.

#### AE design to identify the dimensionality of a suitable latent space.

Given a set of OD measurements ($\mathbf{y}$), we develop an AE to identify
the compressed latent space representation ($\mathbf{z}$) of plasma state. We
use $\mathbf{y^{\prime}}$ to denote the reconstucted input from the AE. In
other words, the encoder maps ${\mathbf{y}}\to{\mathbf{z}}$, and the decoder
maps ${\mathbf{z}}\to{\mathbf{y^{\prime}}}$. In our case, $\mathbf{y}$ (and
hence also $\mathbf{y^{\prime}}$) has 34 values — it comprises four scalar
values ($J_{\parallel}^{\text{sat}}$, $T_{\text{e,u}}$, $N_{\text{t}}$, and
$T_{\text{e,t}}$) and a 30-element 1D array ($P_{\text{rad}}$). However, the
different values (scalars vs. the array) exhibit different ranges of values.
Therefore, we utilize a scaling factor $\alpha_{i}$ for each of the 34 input
dimensions to prevent overfitting on the 30 values for the array. In
particular, we choose $\alpha_{i}=0.1$ for values corresponding to the array
elements and $\alpha_{i}=1.0$ for the scalars, effectively asking the model to
consider each array element 10$\times$ less important than the scalars.
Formally, our AE is trained to minimize the loss function

$\mathcal{L}({\mathbf{y}},{\mathbf{y^{\prime}}})=\frac{1}{N}\sum_{i}\alpha_{i}|y_{i}-y^{\prime}_{i}|^{n},$
(1)

where $N=34$ is the total number of elements in the input OD data
$\mathbf{y}$. For $n=1$, the loss function is the mean absolute error (MAE) or
$L1$ norm, and for $n=2$, it becomes the mean square error (MSE) or $L2$ norm.

Since our input data is simply a collection of values, our AE design uses a
series of fully-connected neural network layers. Each such layer provides
dense connections between input and output values through a linear
transformation (a matrix multiplication with an additive bias) followed by a
nonlinear activation. Given a chosen input/output dimensionality and an
activation function, the training of the model then learns the appropriate
matrix and bias to perform this transformation. We experimented with different
depths of the neural network (i.e., the number of such fully-connected layers)
and the corresponding dimensionalities. All models were trained using the
Tensorflow (Abadi et al., 2015) and Keras (Chollet et al., 2015) frameworks
with the Adam optimizer (Kingma & Ba, 2014). A dataset of 111,598 samples was
used for developing the AE, with a randomly selected 80% subset of this data
for training, and the remaining 20% for validating the model (i.e., assess its
accuracy).

The final AE model, which provided the best reconstruction accuracy (see Table
2), contains three fully-connected layers that progressively reduce the data
dimensionality as $34\to 18\to 10\to 6$. Each layer is followed by the sigmoid
activation function. The output of the last layer is the desired latent space
representation; in other words, a 6-dimensional latent space is adequate to
compress a 34-element diagnostic data in our case so as to recover the input
with sufficient accuracy.

#### VAE design to identify a suitable latent space.

The AE setup described above is capable of generating compressed encoding of
the input data. However, such AEs provide no control of the distribution of
the data in the learned latent space representation. Typically, this is not a
big issue if the AE is used to encode/decide the data, i.e., use the AE (both
encoder and decoder) as one model. However, it is not suitable for our
application as our goal is to utilize only the decoder and in conjunction with
a forward, surrogate model. Specifically in our case, it is important to
provide smoothness guarantees on the latent space representation, otherwise
even small prediction errors from the forward model could be amplified
significantly by the decoder. To regularize the distribution in the latent
space and to improve decoder’s generative capability, we use a specific type
of AE, called $\beta$-variational autoencoder (VAE) (Higgins et al., 2016).
The key difference between the (vanilla) AE and the VAE is that the VAE forces
the distribution of $\mathbf{z}$ to a given reference distribution —
typically, a multivariate normal distribution — offering some important
mathematical guarantees, such as smoothness. To ensure that ${\mathbf{z}}$
mimics the reference distribution, a regularization term is added to the loss
function to measure how different the two distributions are. Specifically, the
loss function from Equation 1 is modified to

$\mathcal{L}({\mathbf{y}},{\mathbf{y^{\prime}}})=\frac{1}{N}\sum_{i}\alpha_{i}|y_{i}-y^{\prime}_{i}|^{n}\quad+\quad\beta
D_{KL}(f({\mathbf{z}})~{}||~{}f_{0}),$ (2)

where $D_{KL}(f({\mathbf{z}})~{}||~{}f_{0})$ is the Kullback-Leibler
divergence, which measures how the distribution $f({\mathbf{z}})$ differs from
the reference distribution $f_{0}$, and a hyperparameter $\beta$ is used to
balance the decoder reconstruction accuracy (the first term in the loss
function) and the orthogonality of latent space coordinates (the second term).
As is common for VAE, our reference distribution $f_{0}$ is a multivariate
normal distribution, $\mathcal{N}(0,1)$. In many applications, $\beta$ is set
to be larger than unity for better separation of independent latent space
coordinate; in our case, forcing LSR distribution function to match normal
distribution is not necessary. Instead, a small $\beta=10^{-9}$ is found to be
suitable. The VAE used in this study has the same network architecture and
training parameters as the AE described above.

###  4.2 Predicting Simulation Outputs using a Surrogate Model

Given the latent space representation, $\mathbf{z}$, generated by the VAE
described above, we next develop another ML model that can predict
$\mathbf{z}$ based on the physics model inputs, $\mathbf{x}$ (i.e., global
discharge parameters). Due to a small numbers of model inputs (i.e., three
varying parameters in 1D UEDGE model) and outputs (i.e., six-dimensional
latent space), here a simple yet robust multilayer perception (MLP) is chosen
to be our forward model. This MLP has four fully-connected layers with 24
neurons on each layer, uses the rectified linear (ReLu) activation function,
and MSE as the loss function. Adam optimizer with a learning rate $0.001$ is
used for training the model.

##  5 Data-driven model performance

In this section, we discuss the results on the performance and validation of
the ML models developed using the methods discussed in section 4.

###  5.1 Latent Space Identification using Autoencoders

Through the AE described above, we identify the dimensionality of a suitable
latent space, $D_{\mathbf{z}}$, to be 6. Table 2 indicates that with both $L1$
or $L2$ norm, the validation loss (in normalized units) increases as
$D_{\mathbf{z}}$ decreases. The $L2$ norm appears to have a critical turning
point at $D_{\mathbf{z}}=6$; the validation loss increases substantially when
$D_{\mathbf{z}}$ is set to be less than $6$, while for $D_{\mathbf{z}}>6$, the
improvement is less profound.

$D_{\mathbf{z}}$ | 8 | 7 | 6 | 5 | 4
---|---|---|---|---|---
$\mathcal{L}({\mathbf{y}},{\mathbf{y^{\prime}}})^{n=1}$ | $3.5\times 10^{-4}$ | $4.2\times 10^{-4}$ | $4.7\times 10^{-4}$ | $4.8\times 10^{-4}$ | $5.2\times 10^{-4}$
$\mathcal{L}({\mathbf{y}},{\mathbf{y^{\prime}}})^{n=2}$ | $2.5\times 10^{-6}$ | $3.1\times 10^{-6}$ | $3.3\times 10^{-6}$ | $4.5\times 10^{-6}$ | $5.2\times 10^{-6}$
Table 2: To choose a suitable dimensionality, $D_{\mathbf{z}}$, of the latent
space, we study the reconstruction errors posed by differnet dimensionalities.
Given the significant jump in the error when going from 6 to 5,
$D_{\mathbf{z}}$ was found to be the appropriate choice. Figure 7: Performance
of $\beta$-VAE in terms of absolute residual error for (a)
$J_{\parallel}^{\text{sat}}$, (b) $n_{\text{e,t}}$, (c) $T_{\text{e,u}}$, (d)
$T_{\text{e,t}}$, (e) $P_{\text{rad}}^{\text{max}}$, and (f) peak radiation,
or detachment front location. Here, $x$-axis is UEDGE (true) value, and the
residual error is defined as
$\varepsilon=f_{\beta\text{-VAE}}-f_{\text{UEDGE}}$ for the quantity $f$. All
predictions show excellent correlation with the true values, as quantified
with almost-perfect $R^{2}$ scores. The distribution of data points are shown
in adjoining panels to overcome the issue of overplotting in the scatter
plots; these indicate that for most cases, relatively low and unbiased error
is obtained, whereas higher errors are associated with relatively smaller
number of samples.

Next, we study the performance of the VAE using Figure 7 with UEDGE data
$f_{\text{UEDGE}}$ on the horizontal axis and the residual between VAE-
reproduced data and UEDGE data
$\varepsilon=f_{\beta\text{-VAE}}-f_{\text{UEDGE}}$ on the vertical axis. If a
perfect reconstruction can be achieved, all data points would lie exactly on
the horizontal $\varepsilon=0$ lines. In our case, two of the four scalar
diagnostic measurements, the electron density at divertor $n_{\text{e,t}}$ and
the electron upstream temperature $T_{\text{e,u}}$ can be reproduced with
excellent accuracy, whereas the other two scalar measurements, the ion
saturation current, $J_{\parallel}^{\text{sat}}$, and the electron temperature
at divertor, $T_{\text{e,t}}$, also show good quality reconstruction, with a
slight performance degradation at high values of the respective diagnostics.
Nevertheless, the $R^{2}$ scores (or, coefficient of determination that
measures the goodness-of-fit between the two variables
$f=\\{f_{i}\\},g=\\{g_{i}\\}$ defined as
$R^{2}=1-\sum_{i}(f_{i}-g_{i})^{2}/\sum_{i}f_{i}^{2}$ with $i$ being the index
of the variable; hence $R^{2}=1$ indicates perfect correlation while $R^{2}=0$
means no correlation) for these four scalars all exceed $0.997$, evidencing a
close to perfect replication of input data from trained $\beta$-VAE
statistical-wise. To quantify the reconstruction quality of the radiation
profile, $P_{\text{rad}}$, we choose two metrics: peak amplitude and location.
The peak radiation amplitude prediction by the VAE is consistent the UEDGE
data with a slightly larger variance ($R^{2}\simeq 0.994<0.997$). The peak
radiation, or detachment front location appears to be challenging to be
correctly reproduced with $R^{2}_{\text{detached}}\approx 0.7$ and
$R^{2}_{\text{attached}}\approx 0$. However, due to the discrete nature of the
simulation mesh, radiation peak locations are discretized as well so that the
$R^{2}$ score here could be misleading, especially for the attached cases. In
fact, all 55,551 attached cases whose peaking radiation location at the
divertor target are correctly reproduced by $\beta$-VAE but they are
represented by what appears to be “a single” orange dot in Figure 7(f).

Figure 8: Residual plots for (a) $T_{e,t}$ near detachment transition, and (b)
$J_{\parallel}^{\text{sat}}/n_{e,t}$. Detached cases are in blue while
attached cases are in orange. Note here $y$-axis of the distribution plots is
in logarithmic scale

One plasma quantity that we are especially interested in for detachment
prediction is the electron temperature at the divertor target,
$T_{\text{e,t}}$, as this is arguably the most important and direct indicator
of detachment. As shown in Figure 8(a), VAE is able to reproduce
$T_{\text{e,t}}$ near the detachment onset value ($\sim 2.1~{}eV$ for this
study) very well. Figure 8(b) shows that the VAE is able to also accurately
retain the ion exist velocity, $v_{\parallel
i,e}=J_{\parallel}^{\text{sat}}/(en_{\text{e,t}})$, which is an inherent
quantity that was not trained directly. A discrepancy between the VAE’s
predictions and UEDGE results appears for
$J_{\parallel}^{\text{sat}}/(en_{\text{e,t}})>55~{}km/s$, which corresponds to
attached plasma cases with $T_{\text{e,t}}\gg 2.1~{}eV$. However, this type of
discrepancy exists only for less than $2\%$ of the total cases - the small
number of $T_{\text{e,t}}\gg 2.1~{}eV$ samples may cause the relatively large
prediction error for these cases. This result indicates that our trained VAE
captures the correct physics constraints between different plasma variables
(i.e., autoencoder inputs and outputs) for at least the majority of the cases
in our dataset.

Finally, we show that $\beta$-VAE encoded LSRs are indeed apart in latent
space based on plasma states (e.g., detached or attached). Figure 9
illustrates the distinctive clusters of LSRs for detached and attached cases
in latent space using t-SNE (t-distributed stochastic neighbor embedding)
method (Van der Maaten & Hinton, 2008). Similarly, Figure 10 has the
distributions of all six latent coordinates, again separated by
attachment/detachment label. The result shows that despite overlaps between
detached and attached plasmas in some latent variables (e.g., latent variables
1, 4 and 5), the two scenarios are broadly well separated in the other latent
variables (e.g., latent variables 2 and 6), which is a nice feature to have
for the upcoming forward detachment prediction.

Figure 9: t-SNE visualisation of a subset of 10,000 randomly picked LSR
samples in the training data set in which detached (blue) and attached
(orange) cases are clearly separated.

###  5.2 Predicting Simulation Outputs using a Surrogate Model

Here, we study the performance of our surrogate model - an MLP - after 10,000
epochs. In particular, we compare the latent space representation generated by
the VAE, $\mathbf{z}$, and that predicted by the surrogate model,
$\mathbf{z_{f}}$. Figure 10 depicts plots the six coordinates within the two
quantities ($\mathbf{z}$ vs. $\mathbf{z_{f}}$) as individual scatter plots and
shows high degree of correlation. As shown in the figure, the $R^{2}$ values
for all six coordinates exceeds $0.99$, indicating exceptional performance and
that our surrogate model is capable of predicting the latent space
representation with almost full accuracy.

Figure 10: Performance of the MLP model in terms of absolute residual error
for all six latent variables. Here $x$-axis is $\beta$-VAE encoded value,
while the residual is defined as
$\varepsilon=f_{\text{MLP}}-f_{\beta\text{-VAE}}$ for quantity $f$. As with
earlier results, an almost-perfect R2 score indicates exceptional prediction
quality.

###  5.3 Forward prediction

With both $\beta$-VAE and forward model trained, we are now able to predict
diagnostic measurement of a plasma state from the model inputs by combining
the forward model and decoder (Figure 6(c)). To ensure that this model is
evaluated properly, a separated validation data set consists 57,655 cases
generated independently following the similar data generation process
described in Section 3.2.

Figure 11: Performance of combined forward prediction model in terms of
residuals for (a) $J_{\parallel}^{\text{sat}}$, (b) $n_{e,t}$, (c) $T_{e,u}$,
(d) $T_{e,t}$, (e) $P_{\text{rad}}^{\text{max}}$, and (f) peak radiation, or
detachment front location. Here $x$-axis is UEDGE (true) value, while the
residual is defined as $f_{\text{residual}}=f_{ML}-f_{\text{UEDGE}}$ for
quantity $f$. The distribution of $f$ and the $R^{2}$ scores between $f_{ML}$
and $f_{\text{UEDGE}}$ for detached (blue) and attached (orange) cases are
also provided.

The performance of forward detachment prediction model is illustrated in
Figure 11 where all 57,655 cases are evaluated with the machine learning model
then validated with the UEDGE simulation results. As expected, the overall
accuracy degrades marginally comparing to the performance of $\beta$-VAE
(e.g., Figure 7) in this sequential MLP and decoder architecture.

Figure 12 and Table 3 elucidate the accuracy of this combined forward
detachment prediction model, quantified with error statistic analysis. We
examine model accuracy for both detached and attached cases, and find that our
model does perform slightly different for different plasma states. For
attached plasma, it gives better results when predicting certain plasma
quantities, such as ion saturation current and divertor target electron
density; while it is less satisfactory for divertor target electron
temperature prediction. However, no significant quantitative difference is
found except the detachment front location prediction. Out of four scalar
measurements, upstream temperature $T_{e,u}$ is the most accurately predicted
diagnostic quantity with mean and standard derivation of relative error
$\mu=-0.37\%$ and $\sigma=0.83$ for the entire validation data set. The other
three scalar measurements are also well predicted with $|\mu|<3\%$ and
$\sigma<4\%$ for all the relative errors. The peak radiation strength
prediction appears to have a fairly large uncertainty ($\sim 10\%$). This is
possibly due to the shape peak structure (i.e., lack of resolution) for
detached plasmas and boundary effect (i.e., maximum value resides at the last
point) for attached plasmas as shown in Figure 13(b). Fortunately, the more
useful information, the peak radiation location (equivalent to the detachment
front location $L_{DF}$ to some extent) prediction is remarkably accurate.
$L_{DF}$ is correctly predicted to be at the divertor target plate for nearly
all attached plasmas, and has a small uncertainty ($\sigma=0.60~{}cm$) for
detached plasmas. Note that the error distributions of four scalar
measurements appear to follow Gaussian distribution, while for the peak
radiation amplitude and location, the error distribution are no long Gaussian
– the standard derivation $\sigma$ is skewed by rare (i.e., probability $\sim
0.1\%$) extreme cases.

Figure 12: Probability distribution functions (PDFs) of relative error of (a) $J_{\parallel}^{\text{sat}}$, (b) $n_{e,t}$, (c) $T_{e,u}$, (d) $T_{e,t}$, (e) $P_{\text{rad}}^{\text{max}}$, and (f) absolution error of peak radiation, or detachment front location. Detached cases are in blue while attached cases are in orange. Dashed lines are shifted Gaussian distribution function with $\mu$ and $\sigma$ from Table 3 correspondingly. | detached | attached | combined
---|---|---|---
| $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
$\epsilon_{J_{\parallel}^{\text{sat}},rel}(\%)$ | 1.66 | 3.76 | 1.67 | 1.26 | 1.67 | 2.79
$\epsilon_{n_{e,t},rel}(\%)$ | 1.06 | 3.97 | 1.81 | 1.66 | 1.44 | 3.05
$\epsilon_{T_{e,u},rel}(\%)$ | -0.43 | 1.05 | -0.30 | 0.52 | -0.37 | 0.83
$\epsilon_{T_{e,t},rel}(\%)$ | -0.13 | 3.90 | -2.70 | 3.06 | -2.01 | 3.57
$\epsilon_{P_{rad}^{max},rel}(\%)$ | 0.19 | 23.76 | 13.39 | 43.85 | 6.88 | 36.00
$\epsilon_{L_{DF},abs}(cm)$ | 0.07 | 0.60 | 0 | 0 | 0.04 | 0.42
Table 3: Mean ($\mu$) and standard deviation ($\sigma$) of relative or
absolute error of forward prediction model. Figure 13: (a) Electron
temperature at divertor target $T_{e,t}$ prediction vs UEDGE simulation
result, and (b) examples of radiation profile or detachment front prediction.
True (UEDGE) detached cases are in blue while attached cases are in orange.

We are particularly interested in $T_{e,t}$ prediction as this is the primary
indicator of detachment. If we use the same detachment criterion (e.g.,
$T_{e,t}=2.1~{}eV$), there are 35 out of 28,468 detached cases (i.e., blue
triangles at the top left quadrant of Figure 13(a)) mis-predicted as
“attached”. Similarly, 25 out of 29,187 attached cases (orange diamonds at the
bottom right quadrant of Figure 13(a)) are mis-labeled “detached”. Even
counting in the other marginal cases, the mis-classification rate is lower
than $0.2\%$. These mis-classified cases appears to congregate near the
detachment onset point (i.e., all these cases have $T_{e,t}\in(1,4)~{}eV$).
However, no common feature can be identified in terms of the control
parameters ${\mathbf{x}}$.

Besides accuracy, speed is another important metric for a real-world
application. Table 4 summarizes our predictive model’s performance in terms of
the wall-clock time required to predict diagnostic measurements based on
controlled inputs. Even without any optimization, this model takes about
$36~{}ms$ to carry out a prediction for one case - about 10,000 times faster
than 1D UEDGE simulations which normally require a few minutes to find a
converged solution. It is already within the minimum requirement for
detachment control ($\sim 100~{}ms$). The model efficiency also increases when
predicting multiple cases (for $n<100,000$).

number of cases | 1 | 10 | 100 | 1,000 | 10,000 | 100,000
---|---|---|---|---|---|---
speed ($ms$) | 36 | 37 | 42 | 51 | 178 | 2,259
Table 4: Wall-clock time v.s. number of cases.

##  6 Comparison with two-point models

Validation of the forward model for detachment prediction shows that our model
is able to accurately predict UEDGE result with a dramatic speed up, but
perhaps a more meaningful benchmark exercise would be comparing our newly
developed model, as well as the UEDGE model with the widely used detachment
prediction models nowadays, namely the analytical basic two-point model (2PM)
(Stangeby et al., 2000), and the most sophisticated semi-analytical two-point
model formatting (2PMF) (Stangeby, 2018).

###  6.1 Basic two-point model

The basic two-point model (2PM) is derived to evaluate $T_{e,u}$, $T_{e,t}$
and $n_{e,t}$ for given $n_{e,u}$, parallel heat flux $q_{\parallel}$ and flux
tube length $L$ based on particle, pressure and power balance. Because of its
simplicity, basic 2PM has been implemented in the tokamak detachment control
as a priori to estimate the degree of detachment (Eldon et al., 2022). Because
2PM doesn’t account for momentum and power loss between the upstream and
downstream points, we therefore benchmark UEDGE, the forward detachment model
and the 2PM on a case with zero-impurity to remove the impurity radiation.
Figure 14 displays upstream density scan of $J_{\parallel}^{\text{sat}}$ and
$T_{e,t}$ for the three models. Clearly, the data-driven model prediction and
the UEDGE simulation results are in good agreement while the 2PM gives a quite
different result. This one order of magnitude discrepancy between our data-
driven model/UEDGE and 2PM is due to the more comprehensive physics in our
model/UEDGE. For instance, in this comparison the upstream ion temperature is
much higher than the electron temperature in the upstream region in the UEDGE
simulations due to slow ion parallel thermal conduction whereas the electron
and ion temperatures are assumed to be the same in the basic 2PM.
Additionally, in the basic 2PM, the parallel heat flux is assumed to be
conductive using the Spitzer-Härm formula; while in 1D UEDGE simulations, the
parallel heat flux is assumed to be local flux-limited thermal transport.
UEDGE is only able to recover 2PM result by further simplifying simulation
(e.g., using slab configuration and turn off non-ideal terms) and manually
adjusting boundary conditions to match 2PM assumptions, such as $T_{e}=T_{i}$
at upstream.

Figure 14: Upstream density scan of (a) $J_{\parallel}^{\text{sat}}$ and (b)
$T_{e,u}$ for UEDGE (blue), machine learning data-driven model (orange) and
basic two-point model(black dashed) with $P_{\text{inj}}=1.71~{}MW$,
$f_{Z}=0$.

###  6.2 Two-point model formatting (2PMF)

To address the volumetric momentum and power losses, and the magnetic geometry
effects, the basic two-point model has been extended to semi-analytical two-
point model formatting (2PMF) with pre-fitted power and momentum loss
coefficients $f_{\text{cooling}}$, $f_{\text{mom-loss}}$ (Stangeby, 2018).
Benchmark between our data-driven model, UEDGE simulation and 2PMF are also
performed. Two fitting curves (from (Stangeby, 2018)) are proposed for both
power and momentum loss coefficients in 2PMF,

Stangeby fitting formula 1:

$\displaystyle 1-f_{\text{cooling}}$
$\displaystyle=[1-\exp({-T_{e,t}/2.4})]^{1.9}$ (3) $\displaystyle
1-f_{\text{mom-loss}}$ $\displaystyle=[1-\exp({-T_{e,t}/0.8})]^{2.1}$ (4)

Stangeby fitting formula 2:

$\displaystyle 1-f_{\text{cooling}}$
$\displaystyle=0.9[1-\exp({-T_{e,t}/6})]^{1.7}$ (5) $\displaystyle
1-f_{\text{mom-loss}}$ $\displaystyle=1.3[1-\exp({-T_{e,t}/1.8})]^{1.6}$ (6)

where UEDGE observed $T_{e,t}$ is used to estimate $f_{\text{cooling}}$ and
$f_{\text{mom-loss}}$. Figure 15 displays upstream density scan of
$J_{\parallel}^{\text{sat}}$ and $T_{e,t}$ for the three models. The results
from 2PMF compare with UEDGE simulation results better than those from the
basic 2PM. All three models show two common features: (1) the ion saturation
current density rollover and (2) temperature cliff at the onset of detachment.
Once again, our data-driven model prediction matches UEDGE simulation results
very well. Although there is no qualitative disagreement between our data-
driven model/UEDGE and 2PMF, quantitatively the predictions between data-
driven model/UEDGE and 2PMF can be off by an order of magnitude. This is
because 2PMF prediction depends heavily on the choice of fitting curves for
the power and momentum loss coefficients which may simplify the nonlinear
dynamics setting the divertor plasma conditions. As illustrated in Figure 15,
fitting formula 1 does a better overall job than fitting formula 2 for this
test case as $T_{e,t}$ predicted by fitting formula 2 is near an order of
magnitude lower than UEDGE result when divertor plasma is detached
($n_{e,u}>2.4\times 10^{19}~{}m^{-3}$).

Figure 15: Upstream density scan of (a) $J_{\parallel}^{\text{sat}}$ and (b)
$T_{e,u}$ for UEDGE (blue), machine learning data-driven model (orange) and
two-point model formatting (black lines) with $P_{\text{inj}}=2.66~{}MW$,
$f_{Z}=0.02$.

##  7 Model applications

As mentioned in the introduction, the main motivation of developing such a
surrogate model is to enable reliable and fast detachment prediction for
integrated machine design, scenario development and real time plasma control.
Heat exhaust is not a severe issue for current tokamaks due to the overall
limited power output. Therefore, the scrape-off layer (SOL) and divertor
plasma dynamics is less of a concern and most of the effort so far has been
focused on exploring operation scenarios with improved fusion performance,
equilibria, and plasma stability control inside the separatrix, e.g., core-
edge integration. However, for future high fusion power devices such as
reactors, their operation space must also fulfill constraints posed by
divertor’s material and engineering limits. Likewise, fusion burn control
needs to incorporate divertor heat and particle exhaust solutions such as
detachment. In other words, designing and operating future devices require
core-edge-SOL/divertor integration. SOL/divertor modelling suffers from the
accuracy-speed trade-off similar to many other research. Current SOL/divertor
transport codes such as UEDGE and SOLPS, are too slow for these applications
as they are designed for physics investigations; while the basic two-point
model and two-point model formatting are fast but perhaps over-simplify the
problem. Our data-driven model overcomes this accuracy-speed gap with some
room to prioritize one factor over another depending on the application. For
instance, speed is the top criterion for plasma control (either in a simulator
or in the actual plasma control system). Real time or even fast than real time
prediction is required in order to activate actuators in time. On the other
hand, accuracy is likely weighted more than speed for device design and
scenario development.

The accuracy of the proposed data-driven approach relies on the quality and
quantity of the training data set. We would like to remark that tokamak edge
plasma contains very rich physics with many factors have influence on the
detachment onset or threshold, such as 2D/3D effects, divertor plate geometry,
multi-species and/or multi-charge-state impurity, wall condition, and so on.
Because of the 1D flux-tube mesh simplification, 1D UEDGE simulations under-
estimate the detachment onset temperature $T_{e,t}$ and over-estimate the peak
radiation amplitude, our current model no surprisingly picks up these
unfavorable predictions. One would expect that the performance of model
improves in terms of matching real experimental measurement when the training
data is extended to incorporate richer physics in a more realistic experiment
setting, i.e., once trained upon higher quality data sets either from higher
fidelity numerical models such as 2D UEDGE/SOLPS-ITER or experiments. Even
though the underlying methodology will be the same, the architecture of neural
networks likely needs to be modified. Notably, synthetic, or real diagnostic
measurements will have a complicated format with 2D simulations or
experiments. In a realistic tokamak, both Langmuir probe and Thompson
scattering are multi-channel; therefore $J_{\parallel}^{\text{sat}}$,
$T_{e,u}$, $n_{e,t}$ and $T_{e,t}$ become sparse 1D arrays in space (point
data), spectroscopic and bolometer diagnostics provide radiation
power/strength at certain wavelength or over the entire spectrum (e.g., 1D
volume averaged data), visible and IR cameras give image (i.e., 2D data with
projected volume and certain range of wavelength averaged quantity). Moreover,
there are likely some discrepancies between model produced synthetic
diagnostic data and real experimental measurements due to various reasons such
as model simplification, inherent instrumentation noise and so on. Handling
these multi-modal diagnostics together in a consistent manner, as well as
bridging the gap between simulation and experimental data will be addressed in
our future study.

Even though the data-driven detachment prediction model presented here is
based on 1D UEDGE simulations, it could be readily integrated for detachment
control application. At first glance, current model is similar to the
analytical basic two-point model which has recently been implemented in KSTAR
tokamak detachment control system (Eldon et al., 2022) as both models are
meant to predict divertor or downstream plasma state. Comparing to the basic
two-point model, our model features additional detachment front prediction.
Additionally, due to fewer simplifications made in UEDGE model such as $B$
variation, local flux-limited thermal transport, impurity radiation effects,
etc, our model should give more reliable predictions than two-point model in
certain circumstances, resulting improved control performance of detachment
and overall plasma confinement.

##  8 Summary

In this paper, we explore a new physics model-based approach to predict
divertor detachment by leveraging the “latent feature space” concept in
machine learning research. As a proof of concept study, highly efficient 1D
UEDGE model which contains the crucial physics ingredients of detachment are
used to simulate the plasma and neutrals along the open magnetic field lines
in the scrape-off layer (SOL) and to generate our training and validation data
sets. Over 160,000 simulations with three varying UEDGE model inputs
${\mathbf{x}}$ (e.g., different upstream density, injection power and carbon
fraction) are performed to cover the normal DIII-D tokamak operation parameter
region; and five synthetic diagnostic measurements such as upstream
temperature, electron density, temperature and saturation current at divertor
target, as well as radiation profile are collected as the diagnostic set
${\mathbf{y}}$. The latent space as well as the latent space representation of
plasma state ${\mathbf{z}}$ are then identified by compressing ${\mathbf{y}}$
through an autoencoder. Sequentially, a forward surrogate model is trained to
make predictions of ${\mathbf{z}}$ from UEDGE model inputs ${\mathbf{x}}$;
then the trained decoder is used to reconstruct diagnostic measurement
${\mathbf{y}}$ back in configuration space. We find that a 6-dimensional
latent space is good enough to closely yield a match for the true system in
configuration space (i.e., the synthetic diagnostic measurements), and the
forward detachment prediction model
(${\mathbf{x}}\to{\mathbf{z}}\to{\mathbf{y}}$) also produces quite accurate
predictions (relative error on the order of a few percent statistically) with
at least $10^{4}$ speed-up comparing to the UEDGE simulations.

Our pilot study demonstrates that the complicated divertor/SOL plasma state
has a low-dimensional representation in latent space. Therefore, this new
latent space description of plasma state not only can be used to construct a
fast and robust surrogate model for steady-state detachment prediction as
revealed in this paper, but also has the potential to be used for dynamical
control once the critical plasma nonlinear dynamics (in latent space) was
successfully identified.

## Acknowledgement

The authors thank the anonymous reviewers for their comments and suggestions.

## Funding

This work is performed under the auspices of the U.S. Department of Energy (US
DOE) by Lawrence Livermore National Laboratory under Contract DE-
AC52-07NA27344 through the pilot project “A critical step for development of
machine learning surrogate models for tokamak divertor-plasma detachment
control". LLNL-JRNL-836030.

## Declaration of interests

The authors report no conflict of interest.

## Data availability statement

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

## References

  * Abadi et al. (2015) Abadi, Martín, Agarwal, Ashish, Barham, Paul, Brevdo, Eugene, Chen, Zhifeng, Citro, Craig, Corrado, Greg S., Davis, Andy, Dean, Jeffrey, Devin, Matthieu, Ghemawat, Sanjay, Goodfellow, Ian, Harp, Andrew, Irving, Geoffrey, Isard, Michael, Jia, Yangqing, Jozefowicz, Rafal, Kaiser, Lukasz, Kudlur, Manjunath, Levenberg, Josh, Mané, Dandelion, Monga, Rajat, Moore, Sherry, Murray, Derek, Olah, Chris, Schuster, Mike, Shlens, Jonathon, Steiner, Benoit, Sutskever, Ilya, Talwar, Kunal, Tucker, Paul, Vanhoucke, Vincent, Vasudevan, Vijay, Viégas, Fernanda, Vinyals, Oriol, Warden, Pete, Wattenberg, Martin, Wicke, Martin, Yu, Yuan & Zheng, Xiaoqiang 2015 TensorFlow: Large-scale machine learning on heterogeneous systems. Software available from tensorflow.org.
  * Anirudh et al. (2022) Anirudh, Rushil, Archibald, Rick, Asif, M Salman, Becker, Markus M, Benkadda, Sadruddin, Bremer, Peer-Timo, Budé, Rick HS, Chang, CS, Chen, Lei, Churchill, RM & others 2022 2022 review of data-driven plasma science. arXiv preprint arXiv:2205.15832 .
  * Anirudh et al. (2020) Anirudh, Rushil, Thiagarajan, Jayaraman J, Bremer, Peer-Timo & Spears, Brian K 2020 Improved surrogates in inertial confinement fusion with manifold and cycle consistencies. Proceedings of the National Academy of Sciences 117 (18), 9741–9746.
  * Ballinger et al. (2021) Ballinger, SB, Kuang, Adam Q, Umansky, MV, Brunner, D, Canik, John M, Greenwald, M, Lore, Jeremy D, LaBombard, B, Terry, James L, Wigram, M & others 2021 Simulation of the SPARC plasma boundary with the UEDGE code. Nuclear Fusion 61 (8), 086014.
  * van den Berg et al. (2018) van den Berg, J, Abramovic, I, Lopes Cardozo, NJ & Moseev, D 2018 Fast analysis of collective thomson scattering spectra on Wendelstein 7-X. Review of Scientific Instruments 89 (8), 083507.
  * Chollet et al. (2015) Chollet, François & others 2015 Keras. https://keras.io.
  * Degrave et al. (2022) Degrave, Jonas, Felici, Federico, Buchli, Jonas, Neunert, Michael, Tracey, Brendan, Carpanese, Francesco, Ewalds, Timo, Hafner, Roland, Abdolmaleki, Abbas, de Las Casas, Diego & others 2022 Magnetic control of tokamak plasmas through deep reinforcement learning. Nature 602 (7897), 414–419.
  * Dudson et al. (2019) Dudson, Benjamin Daniel, Allen, J, Body, T, Chapman, Ben, Lau, C, Townley, Luke, Moulton, D, Harrison, J & Lipschultz, Bruce 2019 The role of particle, energy and momentum losses in 1D simulations of divertor detachment. Plasma Physics and Controlled Fusion 61 (6), 065008.
  * Eldon et al. (2022) Eldon, David, Anand, Himank, Bak, JG, Barr, Jayson L, Hahn, Sang-Hee, Jeong, JH, Kim, H-S, Lee, HH, Leonard, Anthony W, Sammuli, Brian S & others 2022 Enhancement of detachment control with simplified real-time modelling on the KSTAR tokamak. Plasma Physics and Controlled Fusion 64 (7), 075002.
  * Higgins et al. (2016) Higgins, Irina, Matthey, Loic, Pal, Arka, Burgess, Christopher, Glorot, Xavier, Botvinick, Matthew, Mohamed, Shakir & Lerchner, Alexander 2016 beta-VAE: Learning basic visual concepts with a constrained variational framework. ICLR 2017 Conference .
  * Kates-Harbeck et al. (2019) Kates-Harbeck, Julian, Svyatkovskiy, Alexey & Tang, William 2019 Predicting disruptive instabilities in controlled fusion plasmas through deep learning. Nature 568 (7753), 526–531.
  * Kingma & Ba (2014) Kingma, Diederik P & Ba, Jimmy 2014 Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 .
  * Kuang et al. (2020) Kuang, AQ, Ballinger, S, Brunner, D, Canik, John, Creely, AJ, Gray, Travis, Greenwald, M, Hughes, JW, Irby, J, LaBombard, B & others 2020 Divertor heat flux challenge and mitigation in SPARC. Journal of Plasma Physics 86 (5).
  * Loarte et al. (1998) Loarte, Alberto, Monk, RD, Martin-Solis, JR, Campbell, DJ, Chankin, AV, Clement, S, Davies, SJ, Ehrenberg, J, Erents, SK, Guo, HY & others 1998 Plasma detachment in JET Mark I divertor experiments. Nuclear Fusion 38 (3), 331.
  * Ma et al. (2020) Ma, Chenhao, Zhu, Ben, Xu, Xue-Qiao & Wang, Weixing 2020 Machine learning surrogate models for Landau fluid closure. Physics of Plasmas 27 (4), 042502.
  * Van der Maaten & Hinton (2008) Van der Maaten, Laurens & Hinton, Geoffrey 2008 Visualizing data using t-SNE. Journal of Machine Learning Research 9 (11).
  * Mathews et al. (2021) Mathews, Abhilash, Francisquez, Manaure, Hughes, Jerry W, Hatch, David R, Zhu, Ben & Rogers, Barrett N 2021 Uncovering turbulent plasma dynamics via deep learning from partial observations. Physical Review E 104 (2), 025205.
  * Maulik et al. (2020) Maulik, Romit, Garland, Nathan A, Burby, Joshua W, Tang, Xian-Zhu & Balaprakash, Prasanna 2020 Neural network representability of fully ionized plasma fluid model closures. Physics of Plasmas 27 (7), 072106.
  * Maurizio et al. (2017) Maurizio, R, Elmore, S, Fedorczak, N, Gallo, A, Reimerdes, H, Labit, B, Theiler, C, Tsui, CK, Vijvers, WAJ, Team, TCV & others 2017 Divertor power load studies for attached L-mode single-null plasmas in TCV. Nuclear Fusion 58 (1), 016052.
  * McLean et al. (2015) McLean, AG, Leonard, AW, Makowski, MA, Groth, M, Allen, SL, Boedo, JA, Bray, BD, Briesemeister, AR, Carlstrom, TN, Eldon, D & others 2015 Electron pressure balance in the SOL through the transition to detachment. Journal of Nuclear Materials 463, 533–536.
  * Miller et al. (2021) Miller, Marco Andres, Churchill, Randy Michael, Dener, A, Chang, CS, Munson, T & Hager, R 2021 Encoder–decoder neural network for solving the nonlinear Fokker–Planck–Landau collision operator in XGC. Journal of Plasma Physics 87 (2).
  * Montes et al. (2019) Montes, Kevin Joseph, Rea, Cristina, Granetz, RS, Tinguely, Roy Alexander, Eidietis, N, Meneghini, OM, Chen, DL, Shen, Biao, Xiao, BJ, Erickson, Keith & others 2019 Machine learning for disruption warnings on Alcator C-Mod, DIII-D, and EAST. Nuclear Fusion 59 (9), 096015.
  * Peterson et al. (2022) Peterson, J Luc, Bay, Ben, Koning, Joe, Robinson, Peter, Semler, Jessica, White, Jeremy, Anirudh, Rushil, Athey, Kevin, Bremer, Peer-Timo, Di Natale, Francesco & others 2022 Enabling machine learning-ready HPC ensembles with Merlin. Future Generation Computer Systems 131, 255–268.
  * Porter et al. (1996) Porter, Gary D, Allen, SL, Brown, M, Fenstermacher, ME, Hill, DN, Jong, RA, Leonard, AW, Nilson, D, Rensink, ME, Rognlien, TD & others 1996 Simulation of experimentally achieved DIII-D detached plasmas using the UEDGE code. Physics of Plasmas 3 (5), 1967–1975.
  * Rea et al. (2020) Rea, Christina, Montes, Kevin J, Pau, A, Granetz, RS & Sauter, O 2020 Progress toward interpretable machine learning–based disruption predictors across tokamaks. Fusion Science and Technology 76 (8), 912–924.
  * Rognlien et al. (1994) Rognlien, TD, Brown, PN, Campbell, RB, Kaiser, TB, Knoll, DA, McHugh, PR, Porter, GD, Rensink, ME & Smith, GR 1994 2-D fluid transport simulations of gaseous/radiative divertors. Contributions to Plasma Physics 34 (2-3), 362–367.
  * Rognlien & Rensink (2002) Rognlien, TD & Rensink, ME 2002 Edge-plasma models and characteristics for magnetic fusion energy devices. Fusion Engineering and Design 60 (4), 497–514.
  * Samuell et al. (2021) Samuell, CM, Mclean, AG, Johnson, CA, Glass, F & Jaervinen, AE 2021 Measuring the electron temperature and identifying plasma detachment using machine learning and spectroscopy. Review of Scientific Instruments 92 (4), 043520.
  * Stangeby (2018) Stangeby, PC 2018 Basic physical processes and reduced models for plasma detachment. Plasma Physics and Controlled Fusion 60 (4), 044022.
  * Stangeby et al. (2000) Stangeby, Peter C & others 2000 The plasma boundary of magnetic fusion devices, , vol. 224. Institute of Physics Pub. Philadelphia, Pennsylvania.
  * Wang et al. (2020) Wang, Libo, Xu, XQ, Zhu, Ben, Ma, Chenhao & Lei, Yi-an 2020 Deep learning surrogate model for kinetic Landau-fluid closure with collision. AIP Advances 10 (7), 075108\.
  * Wiesen et al. (2015) Wiesen, Sven, Reiter, Detlev, Kotov, Vladislav, Baelmans, Martine, Dekeyser, Wouter, Kukushkin, AS, Lisgo, SW, Pitts, RA, Rozhansky, Vladimir, Saibene, Gabriella & others 2015 The new SOLPS-ITER code package. Journal of Nuclear Materials 463, 480–484.
  * Wising et al. (1996) Wising, F, Knoll, DA, Krasheninnikov, SI, Rognlien, TD & Sigmar, DJ 1996 Simulation of detachment in ITER-geometry using the UEDGE code and a fluid neutral model. Contributions to Plasma Physics 36 (2-3), 309–313.
  * Wising et al. (1997) Wising, F, Krasheninnikov, SI, Sigmar, DJ, Knoll, DA, Rognlien, TD, LaBombard, B, Lipschultz, B & McCracken, G 1997 Simulation of plasma flux detachment in Alcator C-Mod and ITER. Journal of Nuclear Materials 241, 273–277.
  * Zhao et al. (2022) Zhao, Menglong, Rognlien, Tom, Zhu, Ben, Meyer, William, Xu, Xueqiao, Dudson, Ben & Li, Nami 2022 Sub-supersonic transitions in the tokamak scrape-off layer and bifurcation solutions. Being prepared for Physics of Plasmas .