††thanks<EMAIL_ADDRESS>
# Predicted critical state based on invariance of the Lyapunov exponent in
dual spaces

Tong Liu  Department of Applied Physics, School of Science, Nanjing University
of Posts and Telecommunications, Nanjing 210003, China  Xu Xia  Academy of
Mathematics and System Sciences, Chinese Academy of Sciences, Beijing 100190,
China

(September 3, 2024)

###### Abstract

The critical state in disordered systems, a fascinating and subtle eigenstate,
has attracted a lot of research interest. However, the nature of the critical
state is difficult to describe quantitatively, and in general it cannot
predict a system that host the critical state. In this work, we propose an
explicit criterion that Lyapunov exponent of the critical state should be 0
simultaneously in dual spaces, namely Lyapunov exponent remains invariant
under Fourier transform. With this criterion, we exactly predict a system
hosting a large number of critical states. Then, we perform numerical
verification of the theoretical prediction, and display the self-similarity of
the critical state. Finally, we conjecture that there exists some kind of
connection between the invariance of the Lyapunov exponent and conformal
invariance, which can promote the research of critical phenomena.

##  I Introduction

Since the phenomenon of Anderson localization Anderson  was proposed, the
quantum disordered system has attracted extensive research enthusiasm. When
impurities are randomly added to an ideal conductor, Bloch state in the
conductor will change into localized state with the increase of the average
concentration of impurities. Heretofore many landmark achievements have been
obtained scaling1 ; scaling2 . For a three-dimensional tight-binding model
with short-range hopping and random on-site potential, there exists a critical
energy separating extended state and localized state, which is dubbed mobility
edge Mott . And many rigorous mathematical methods were developed to
quantitatively analyze Anderson localization phenomenon A50 . An example is to
use the transfer matrix method to solve the Lyapunov exponent of one-
dimensional Anderson model MA , as one-dimensional problem is more likely to
be solved analytically.

Transfer matrix method focuses on the recurrence relation of the amplitude of
wave functions, and is to calculate the average exponential divergence rate
Soukoulis . Mathematically, let Lyapunov exponent $\gamma$ be rewritten as the
form of $\exp(-\gamma)$, the latter indicates the divergence rate between
adjacent lattice points. In position space, for an extended state, the
amplitude between neighboring lattice points should be equal under the
thermodynamic limit, so $\gamma$ should be 0; whereas for a localized state,
the amplitude between neighboring lattice points should be exponentially
decaying, so $\gamma$ should be greater than 0. Consequently, we can utilize
the Lyapunov exponent to explicitly distinguish the state, namely, the
extended state corresponds to $\gamma=0$ and localized state corresponds to
$\gamma>0$.

However, in the quantum disordered system, there is a class of rare and
important states, namely the critical state, the rarity means that it
generally only appears at the phase transition point Sokoloff , such as the
eigenstate at the mobility edge. The critical state is neither an extended
state nor a localized state, its prominent feature is possessing the self-
similar structure A50 . In the language of multifractal theory fractal , the
minimum scaling index of the extended state should tend to 1, that of the
localized state should tend to 0, and that of the critical state should be
greater than 0 but less than 1. However, Multifractal theory requires that the
system should be numerically diagonalized firstly, then using numerically
obtained eigenstates to calculate the minimum scaling index. This method
cannot predict the nature of the eigenstate before numerically diagonalizing
the Hamitonian of a system.

In contrast, calculating the Lyapunov exponent does not need detailed
information about eigenstates, and it can directly be solved by giving the
Hamiltonian Avila1 . As long as Lyapunov exponent of the system being
obtained, we can predict the properties of the eigenstates without numerical
diagonalization. For the famous Aubry-André model AA , Lyapunov exponent is
$\gamma=\max\\{0,\ln(V/2)\\}$ in position space, and
$\gamma_{m}=\max\\{0,\ln(2/V)\\}$ in momentum space, where $V$ is the strength
of quasiperiodic potential. We can accurately predict when $0<V<2$,
$\gamma=0$, the eigenstate in position space is extended; while $V>2$,
$\gamma>0$, the eigenstate is localized. However, a confusion emerges, when
$V=2$ (the phase transition point), $\gamma$ is also equal to 0, at this point
the eigenstates are all critical states. Consequently, the Lyapunov exponent
cannot distinguish the extended state from the critical state, the $\gamma=0$
method is invalid.

Therefore, numerous research efforts are devoted to accurately describing the
properties of critical states. Reference Goncalves  develops a
renormalization-group theory to describe the localization properties of
quasiperiodic systems. This theory can be used to get exact or approximate
analytical expressions for the phase boundaries of extended, localized and
critical phases for specific models. Reference GM  proposes that the coupling
between localized and extended states in their overlapped spectrum can provide
a general recipe to construct critical states for two-chain models. Reference
LXJ  points out that there are two approaches to involve the existence of
critical states, one is involving unbounded potential LT  and the other is
involving zeros of hopping terms in Hamiltonian. However, an explicit
criterion for characterizing critical states is still lacking.

To fill this theoretical gap, we give a quantitative criterion to predict the
critical state in general sense. Let’s re-examine Lyapunov exponent in both
position space and momentum space, for extended state, we have $\gamma=0$ and
$\gamma_{m}>0$; for localized state, we have $\gamma>0$ and $\gamma_{m}=0$. In
order to distinguish from the above two states while maintaining the
approaching 0 property of Lyapunov exponent, a reasonable conjecture is that
the critical state corresponds to the invariance of Lyapunov exponent in dual
spaces, namely its distinguishing feature is satisfying the condition

$\gamma=\gamma_{m}=0.$ (1)

To verify the validity of this conjecture, we introduce an exactly solvable
model hosting critical states in a wide range of parameters, and exactly
predict the interval of the existence of critical states through
$\gamma=\gamma_{m}=0$. Hence we demonstrate that this conjecture is applicable
to the various critical states.

##  II Model and prediction

The difference equation of the model that we consider is written as,

$\psi_{n+1}+\psi_{n-1}+Vi\tan(2\pi\alpha n+\theta)\psi_{n}=E\psi_{n},$ (2)

where $V$ is the complex potential strength, $E$ is the eigenvalue of systems,
and $\psi_{n}$ is the amplitude of wave function at the $n$th lattice. We
choose to unitize the nearest-neighbor hopping amplitude, a typical choice for
irrational parameter is $\alpha=(\sqrt{5}-1)/2$, and $\theta$ is the phase
factor.

With the Hamiltonian of the system [Eq. (2)], we can determine Lyapunov
exponent $\gamma$ in position space, which is calculated by taking the product
of the transfer matrix $T(\theta)$, namely
$\gamma=\lim_{n\rightarrow\infty}\ln||T_{n}(\theta)||/n.$ Utilizing Avila’s
global theory Avila , we can get the explicit expression of $\gamma$, see
detailed calculation in Supplemental Material Supplemental ,

$\displaystyle\gamma(E)=\max\\{$
$\displaystyle\operatorname{arcosh}\frac{\left|E+V+2\right|+\left|E+V-2\right|}{4},$
(3)
$\displaystyle\operatorname{arcosh}\frac{\left|E-V+2\right|+\left|E-V-2\right|}{4}\\}.$

From Eq. (3), we can extract the allowed energies of the system. Firstly, we
make $\gamma(E)=0$, and obtain that $E$ is within the region $[V-2,2-V]$.
Secondly, we make $\gamma(E)>0$, and obtain that $E$ is within the region
$\\{iy\leavevmode\nobreak\ |\leavevmode\nobreak\ y\in\mathbb{R}^{*}\\}$
($V\leq 2$) or $\\{iy\leavevmode\nobreak\ |\leavevmode\nobreak\
y\in\mathbb{R}\\}$ ($V>2$), which means that when the eigenstate is a
localized state in the position space, the eigenvalue is a pure imaginary
number. However, when $\gamma=0$, with regard to the eigenvalue within
$[V-2,2-V]$, it cannot be concluded that the corresponding eigenstate is an
extended state or a critical state.

The next step is to calculate Lyapunov exponent $\gamma_{m}$ in momentum
space. Firstly, we introduce the Fourier transformation

$f_{k}=\frac{1}{\sqrt{L}}\sum_{n=1}^{L}e^{i2\pi\alpha kn}\psi_{n},$

thus the dual equation of Eq. (2) in momentum space is written as

${f_{k+1}}=\frac{-2\cos[2\pi(k-1)\alpha]+V+E}{2\cos[2\pi(k+1)\alpha]+V-E}f_{k-1}.$
(4)

From Eq. (4), an initial wave function solution can be written as

$f_{k}\propto\left\\{\begin{array}[]{cl}0,&k=0,\pm 2,\pm 4\cdots,\\\
0,&k=2j+1<k_{0},\\\ 1,&k=2j+1=k_{0},\\\
\frac{g^{(1)}_{k-2}}{g^{(2)}_{k}}f_{k-2},&k=2j+1>k_{0},\end{array}\right.$

then Lyapunov exponent $\gamma_{m}$ can be obtained by Sarnark’s method
Supplemental ; Sarnak ,

$\gamma_{m}(E)=\lim_{k\rightarrow\infty}\frac{1}{k-k_{0}}\ln\left|\frac{f_{k}}{f_{k_{0}}}\right|=\frac{1}{2\pi}\int_{0}^{2\pi}\ln
g^{(1)}-\ln g^{(2)}d\theta,$ (5)

where $g^{(1)}=|-2\cos(2\pi\theta)+V+E|,g^{(2)}=|2\cos(2\pi\theta)+V-E|$.

Recalling the calculation process of Lyapunov exponent $\gamma$ in position
space, when $\gamma>0$, namely eigenvalues of the system are pure imaginary
numbers, we have the identity

$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}\ln|-2\cos(2\pi\theta)+V+iy|d\theta$
(6)
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\ln|2\cos(2\pi\theta)+V-iy|d\theta,$

hence Lyapunov exponent in momentum space $\gamma_{m}=0$. $\gamma>0$ and
$\gamma_{m}=0$ indicate the associated state is the extended state in momentum
space, which is exactly corresponding to the localized state in position
space.

When $\gamma=0$, namely eigenvalues are real numbers and within the interval
$[V-2,2-V]$, we have

$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}\ln|g^{(1)}|d\theta=\left\\{\begin{array}[]{cl}0,&g^{(1)}=0,\\\
\ln\left|\frac{|E+V|+\sqrt{(E+V)^{2}-4}}{2}\right|,&g^{(1)}\neq
0,\end{array}\right.$ (7)

and

$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}\ln|g^{(2)}|d\theta=\left\\{\begin{array}[]{cl}0,&g^{(2)}=0,\\\
\ln\left|\frac{|E-V|+\sqrt{(E-V)^{2}-4}}{2}\right|,&g^{(2)}\neq
0.\end{array}\right.$ (8)

According to Eq. (7) and Eq. (8), we obtain the identity

$\frac{1}{2\pi}\int_{0}^{2\pi}\ln|g^{(1)}|d\theta=\frac{1}{2\pi}\int_{0}^{2\pi}\ln|g^{(2)}|d\theta,$

see detailed calculation in Supplemental Material Supplemental . Hence, when
$V-2\leq E\leq 2-V$, $\gamma_{m}$ is also equal to 0, which is very similar to
the result that the eigenvalues are pure imaginary numbers. The difference is
that the latter have $\gamma>0$ in position space and $\gamma_{m}=0$ in
momentum space; whereas the eigenstates of the former have
$\gamma=\gamma_{m}=0$.

Based on the above achievements and the conjecture Eq. (1), we can predict
that different from the common models, where the critical state only exists at
the phase transition point, for the model of Eq. (2), critical states exist in
a wide range of the parameters $0<V\leq 2$, due to $V-2\leq E\leq 2-V$
indicating $0<V\leq 2$. Consequently, we realized the prediction of a system
hosting a large number of critical states through giving its Hamiltonian.

##  III Numerical verification and Self-similarity

Figure 1: (Color online) The black curve represents wave function of $E=-1$
obtained from Eq. (2) with the parameter $V=1$. The red curves represent three
wave function peaks after magnifying. Starting with the smallest peak, the
scaled multiples are 7, 4 and 2.2857 in turn. It clearly shows that the scaled
three smaller peaks are very similar to the largest peak. The total number of
sites is set to be $L=6765$.
Figure 2: (Color online) $\beta_{min}$ as a function of the inverse Fibonacci
index $1/m$ for various eigenvalues. In the large size limit, the brown "x"
markers tend to 0.67 and correspond to the critical state with $E=-1$
($V=0.5$), the brown "o" markers tend to 0.58 and correspond to the extended
state with $E=-1$ ($V=1$), the red square markers tend to 0 and correspond to
the localized state with $E=i$ ($V=3$).

To support the analytical results given in the previous section, we perform
numerical calculations and analyze the critical nature of the eigenstate. We
numerically diagonalize the Hamiltonian (2) in a large size, and get the
eigenvalue and the associated eigenstates. To verify the critical state, we
should illustrate the distinctive feature of the corresponding state, namely
self-similarity.

Mathematically, self-similarity is a typical property of the fractal fractal1
; fractal2 , a self-similar construction is exactly or approximately similar
to a part of itself. Many physical objects in nature, such as capillary
distribution and leaf veins, are statistically self-similar: parts of them
show the same statistical properties as the whole. An equivalent description
of self-similarity is scale invariant fractal3 ; fractal4 , where there is a
smaller part that similar to the proximate larger part at the certain
amplitude. Thus, we can perform the contraction-expansion variations on the
part of physical quantities of the system. If the variation amplitude meets a
certain value, it means that the physical quantity is scale invariant, that
is, the system has a self-similarity structure.

Therefore, we directly numerically diagonalize Eq. (2), plot the obtained
eigenstate, and examine detailed structure of wave function to determine
whether there exists self-similarity. As shown in Fig. 1, the black curve is
the original eigenstate of the system, it has four wave-function peaks of
different amplitudes. This state is not a localized state, but it does not
look like an extended state, because the peaks of the extended state should be
equal. When we magnify three smaller peaks by 7, 4 and 2.2857 times
successively, namely the red curves, they become very similar to the largest
black peak. More importantly, we found that $7/4\approx 4/2.2857\approx 1.75$,
which indicates that three smaller peaks can be expanded to the largest one by
a certain multiple $1.75$. This clearly shows that wave function of the system
has scale invariance and self similarity, hence this state is definitely a
critical state. We also calculate the eigenstates corresponding to different
energy levels and different sizes, and all numerical results are critical
states as expected, see detailed calculation in Supplemental Material
Supplemental . For the critical state of known models known1 ; known2 , the
validity of the criterion can also be conveniently verified.

In addition to visually displaying the wave function, we also calculate the
minimum scaling index of the critical state according the multifractal theory
fractal . For giving the wave function $\psi_{n}$, a scaling index $\beta_{n}$
can be extracted from the $n$th on-site probability
$P_{n}=|\psi_{n}|^{2}\sim(1/F_{m})^{\beta_{n}}$, where $F_{m}$ is the $m$th
Fibonacci number. The multifractal theorem states, when the wave functions are
extended, the maximum of $P_{n}$ scales as $\max(P_{n})\sim(1/F_{m})^{1}$,
i.e., $\beta_{min}=\min(\beta_{n})=1$. On the other hand, when the wave
functions are localized, $P_{n}$ concentrates at the individual site and tends
to zero at the other sites, yielding $\max(P_{n})\sim(1/F_{m})^{0}$ and
$\beta_{min}=\min(\beta_{n})=0$. With regard to the critical state, the
corresponding $\beta_{min}$ is located within the interval
$\left(0,\leavevmode\nobreak\ 1\right)$, and can be utilized to distinguish
extended and critical states. In order to reduce finite-size effects, we
examines the trend of $\beta_{min}$ under the limit of large size.

As shown in Fig. 2, $\beta_{min}$ is plotted as a function of the inverse
Fibonacci index $1/m$, when $m\rightarrow 0$, the system size
$L\rightarrow\infty$. It clearly shows that $\beta_{min}$ is between $0$ and
$1$ in the large $L$ limit for the eigenvalues $E=-1$ ($V=0.5$) and $E=-1$
($V=1$), hence the corresponding state is critical. While for the eigenenergy
$E=i$ ($V=3$), $\beta_{min}$ asymptotically tends to 0 in the large $L$ limit,
indicating that the corresponding state is localized. We have also checked
other combinations of parameters and get the same results as expected. Thus,
above numerical results are in excellent agreement with the analytical
results.

##  IV Conclusion

In this work, we propose a conjecture to predict and identify the existence of
the critical state in a system, that is, Lyapunov exponent of the eigenstate
should be 0 in both position space and momentum space. To illustrate this
criterion, we introduce an exactly solvable model, predict and verify the
existence of a large number of critical states in the wide range of the
potential strength. This demonstrates that $\gamma=\gamma_{m}=0$ is not
limited to the critical state at the phase transition point, but is applicable
to various types of critical states. Our findings provide an explicit
quantitative description of the characteristics of the critical state, and
have a positive significance for exploring the critical state of unknown
systems.

Furthermore, it is noteworthy that the prominent feature of critical states is
the scale invariance of wave functions, which indicates that the invariance of
Lyapunov exponent under Fourier transform has a closed relation with conformal
invariance. This is likely to be a new application scenario of conformal
invariance theory, which can describe properties of quantum disordered systems
near the critical point or within the critical interval.

###### Acknowledgements.

T. L. thanks Ming Gong for beneficial comments. This work was supported by the
Natural Science Foundation of Jiangsu Province (Grant No. BK20200737), NUPTSF
(Grants No. NY220090 and No. NY220208), the Innovation Research Project of
Jiangsu Province (Grant No. JSSCBS20210521), and China Postdoctoral Science
Foundation (Grant No. 2022M721693).

## References

  * (1) P. W. Anderson, Absence of diffusion incertain random lattices, Phys. Rev. 109, 1492 (1958).
  * (2) E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Scaling theory of localization: Absence of quantum diffusion in two dimensions, Phys. Rev. Lett. 42, 673 (1979).
  * (3) L. Fleishman and D. C. Licciardello, Fluctuations and localization in one dimension, J. Phys. C 10, L125 (1977).
  * (4) N. Mott, The mobility edge since 1967, J. Phys. C 20, 3075 (1987).
  * (5) A. Lagendijk, B. van Tiggelen, and D. S. Wiersma, Fifty years of Anderson localization. Phys. Today 62, 24 (2009).
  * (6) T. Brandes and S. Kettemann, The Anderson Transition and its Ramifications — Localisation, Quantum Interference, and Interactions. (Springer, Berlin, 2003).
  * (7) C. M. Soukoulis and E. N. Economou, Localization in One-Dimensional Lattices in the Presence of Incommensurate Potentials, Phys. Rev. Lett. 48, 1043 (1982).
  * (8) J. B. Sokoloff, Unusual band structure, wave functions and electrical conductance in crystals with incommensurate periodic potentials, Phys. Rep. 126, 189 (1985).
  * (9) H. Hiramoto and M. Kohmoto, Scaling analysis of quasiperiodic systems: Generalized Harper model, Phys. Rev. B 40, 8225 (1989).
  * (10) A. Avila, J. You , Q. Zhou, Sharp phase transitions for the almost Mathieu operator, Duke. Math. J. 14, 166 (2017).
  * (11) S. Aubry and G. André, Analyticity breaking and Anderson localization in incommensurate lattices, Ann. Israel Phys. Soc. 3(133), 18 (1980).
  * (12) M. Gonçalves, B. Amorim, E. V. Castro, P. Ribeiro, Renormalization-Group Theory of 1D quasiperiodic lattice models with commensurate approximants, arXiv:2206.13549v2.
  * (13) X. Lin, X. Chen, G. Guo, M. Gong, General approach to tunable critical phases with two coupled chains, arXiv:2209.03060v1.
  * (14) X. Zhou, Y. Wang, T. J. Poon, Q. Zhou, X. Liu, Exact new mobility edges between critical and localized states, arXiv:2212.14285v2.
  * (15) T. Liu, X. Xia, S. Longhi, and L. Sanchez-Palencia, Anmalous mobility edges in one-dimensional quasiperiodic models, SciPost Phys. 12, 027 (2022).
  * (16) A. Avila, Global theory of one-frequency Schrödinger operators, Acta. Math. 1, 215, (2015).
  * (17) See Supplemental Material for details of (i) Derivation of Lyapunov exponent $\gamma$ in position space, (ii) Derivation of Lyapunov exponent $\gamma_{m}$ in momentum space, (iii) More numerical verification.
  * (18) P. Sarnak, Spectral behavior of quasi periodic potentials, Commun. Math. Phys. 84, 377 (1982).
  * (19) K. Amin, R. Nagarajan, R. Pandit, and A. Bid, Multifractal Conductance Fluctuations in High-Mobility Graphene in the Integer Quantum Hall Regime, Phys. Rev. Lett. 129, 186802 (2022).
  * (20) X. Deng, S. Ray, S. Sinha, G. V. Shlyapnikov, and L. Santos, One-Dimensional Quasicrystals with Power-Law Hopping, Phys. Rev. Lett. 123, 025301 (2019).
  * (21) H. Yao, A. Khoudli, L. Bresque, and L. Sanchez-Palencia, Critical Behavior and Fractality in Shallow One-Dimensional Quasiperiodic Potentials, Phys. Rev. Lett. 123, 070405 (2019).
  * (22) A. Wardak and P. Gong, Extended Anderson Criticality in Heavy-Tailed Neural Networks, Phys. Rev. Lett. 129 048103 (2022).
  * (23) S. Longhi, Metal-insulator phase transition in a non-Hermitian Aubry-André-Harper Model, Phys. Rev. B 100, 125157 (2019).
  * (24) F. Liu, S. Ghosh, and Y. D. Chong, Localization and adiabatic pumping in a generalized Aubry-André-Harper model, Phys. Rev. B 91, 014108 (2015).

Supplementary Material

This Supplemental Material provides additional information for the main text.
In Sec. S-1, we provide the details of calculating Lyapunov exponent in
position space. In Sec. S-2, we provide the details of calculating Lyapunov
exponent in momentum space. Finally, we give more numerical verification of
critical state in Sec. S-3.

##  S-1 Derivation of Lyapunov exponent $\gamma$ in position space

The Lyapunov exponent $\gamma$ can be calculated by taking the product of the
transfer matrix $T(\theta)$, namely multiplying the transfer matrix $n$ times
consecutively, which is written as

$T_{n}(\theta)=\prod_{l=0}^{n-1}T(2\pi\alpha
l+\theta)=\prod_{l=0}^{n-1}\left(\begin{array}[]{cc}E-Vi\tan(2\pi\alpha
l+\theta)&-1\\\ 1&0\end{array}\right),$

then Lyapunov exponent is $\ln||T_{n}(\theta)||/n$ as $n$ tends to the
infinite in the thermodynamic limit.

The method we use here to calculate Lyapunov exponent is the complexified
phase approach, specifically by continuing the imaginary part of the phase
$\epsilon$, we focus on the new Lyapunov exponent, that is

$T_{n}(\theta+i\epsilon)=\prod_{l=0}^{n-1}T(2\pi\alpha l+\theta+i\epsilon),$

correspondingly, we get $\gamma(\epsilon)$ is
$\lim_{n\rightarrow\infty}\ln||T_{n}(\theta+i\epsilon)||/n.$

Relying on Avila’s global theory Avila2015 , if we can obtain Lyapunov
exponent $\gamma(\epsilon)$ when $\epsilon$ is sufficiently large, then we can
trace back to the specific Lyapunov exponent $\gamma(0)$ when $\epsilon=0$,
namely the original Lyapunov exponent $\gamma$ in position space.

Firstly, rewriting the transfer matrix

$\displaystyle T(\theta)$
$\displaystyle=\left(\begin{array}[]{cc}E-Vi\tan(2\pi\alpha l+\theta)&-1\\\
1&0\end{array}\right)$ (S1) $\displaystyle=\sec(2\pi\alpha l+\theta)B(\theta)$

where

$\displaystyle B(\theta)=\left(\begin{array}[]{cc}E\cos(2\pi\alpha
l+\theta)-Vi\sin(2\pi\alpha l+\theta)&-\cos(2\pi\alpha l+\theta)\\\
\cos(2\pi\alpha l+\theta)&0\end{array}\right)$ (S2)

then, $\gamma(\epsilon)$ can be expressed as

$\displaystyle\gamma_{\epsilon}(E)=$
$\displaystyle\lim_{n\rightarrow\infty}\ln||T_{n}(\theta+i\epsilon)||/n$ (S3)
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}\int\ln|B_{n}(\theta+i\epsilon)|d\theta+\int\ln|\sec(\theta+i\epsilon)|d\theta.$
$\displaystyle=$ $\displaystyle\gamma_{\epsilon}^{B}(E)+\ln(2)-2\pi\epsilon,$

where

$\displaystyle\gamma_{\epsilon}^{B}(E)=\lim_{n\rightarrow\infty}\frac{1}{n}\int\ln\left|B_{n}(\theta+i\varepsilon)\right|d\theta,$
$\displaystyle B_{n}(\theta+ic)=\prod_{l=0}^{n-1}B(2\pi\alpha
l+\theta+i\epsilon).$

When $\epsilon$ tends to $+\infty$, a direct calculating result of
$B(\theta+i\epsilon)$ is

$B(\theta+i\epsilon)=-\frac{1}{2}e^{2\pi\epsilon}e^{i2\pi(\theta+\alpha)}\left(\begin{array}[]{cc}E-V&-1\\\
\\\ 1&0\\\ \\\ \end{array}\right)+o(1).$ (S4)

Thus we get
${\gamma^{+B}}_{\epsilon}(E)=2\pi\epsilon+\ln|\frac{\sqrt{(E-V)^{2}\pm
4}}{2}|-\ln(2)+o(1)$. As a function of $\epsilon$, $\gamma_{\epsilon}^{B}(E)$
is a convex, piecewise linear function whose slope is an integer multiplied by
$2\pi$, hence it is concluded that when $\epsilon$ tends to infinity, we
obtain${\gamma^{+B}}_{\epsilon}(E)=2\pi\epsilon+\ln|\frac{\sqrt{(E-V)^{2}\pm
4}}{2}|-\ln(2)$. And according to the equation(S3), it leads that when
$\epsilon$ is the very large positive number,
${\gamma^{+}}_{\epsilon}(E)={\gamma^{+B}}_{\epsilon}(E)+\ln(2)-2\pi\epsilon=\ln|\frac{\sqrt{(E-V)^{2}\pm
4}}{2}|.$

When $\epsilon$ tends to $-\infty$, a direct calculating result of
$B(\theta+i\epsilon)$ is

$B(\theta+i\epsilon)=-\frac{1}{2}e^{2\pi\epsilon}e^{i2\pi(\theta+\alpha)}\left(\begin{array}[]{cc}E+V&-1\\\
\\\ 1&0\\\ \\\ \end{array}\right)+o(1).$ (S5)

Thus we get
${\gamma^{-B}}_{\epsilon}(E)=2\pi\epsilon+\ln|\frac{\sqrt{(E+V)^{2}\pm
4}}{2}|-\ln(2)+o(1)$. As a function of $\epsilon$, $\gamma_{\epsilon}^{B}(E)$
is a convex, piecewise linear function whose slope is an integer multiplied by
$2\pi$, hence it is concluded that when $\epsilon$ tends to infinity, we
obtain ${\gamma^{-B}}_{\epsilon}(E)=2\pi\epsilon+\ln|\frac{\sqrt{(E+V)^{2}\pm
4}}{2}|-\ln(2)$. And according to equation(S3), it leads that when $\epsilon$
is the very large negative number,
${\gamma^{-}}_{\epsilon}(E)={\gamma^{-B}}_{\epsilon}(E)+\ln(2)-2\pi\epsilon=\ln|\frac{\sqrt{(E+V)^{2}\pm
4}}{2}|.$

Since $\gamma_{\epsilon}$ is a convex function in two semilinear $(0,+\infty)$
and $(-\infty,0)$, it is linear in the cross section and the slope is an
integer multiplied by the $2\pi$ integer, Lyapunov exponent is

$\gamma_{\varepsilon}(E)=\left\\{\begin{array}[]{cc}\gamma_{\varepsilon}^{+}(E)&\varepsilon>0,\\\
\gamma_{\varepsilon}^{+}(E)+2\varepsilon&\frac{\gamma_{\varepsilon}^{-}(E)-\gamma_{\varepsilon}^{+}(E)}{2}<\varepsilon<0,\\\
\gamma_{\varepsilon}^{-}(E)&\varepsilon<\frac{\gamma_{\varepsilon}^{-}(E)-\gamma_{\varepsilon}^{+}(E)}{2},\end{array}\right.$

the relationship between the left and right limit conditions is
$\gamma_{\varepsilon}^{+}(E)>\gamma_{\varepsilon}^{-}(E)$ for any given value
of $\varepsilon$.

Similarly, if the large and small relationship between the left and right
limits is $\gamma_{\varepsilon}^{+}(E)<\gamma_{\varepsilon}^{-}(E)$, Lyapunov
exponent is

$\text{ then
}\gamma_{\varepsilon}(E)=\left\\{\begin{array}[]{cc}\gamma_{\varepsilon}^{-}(E)&\varepsilon<0,\\\
\gamma_{\varepsilon}^{-}(E)-2\varepsilon&0<\varepsilon<\frac{\gamma_{\varepsilon}^{-}(E)-\gamma_{\varepsilon}^{+}(E)}{2},\\\
\gamma_{\varepsilon}^{+}(E)&\varepsilon>\frac{\gamma_{\varepsilon}^{-}(E)-\gamma_{\varepsilon}^{+}(E)}{2}.\end{array}\right.$

Summarizing the above conclusions, Lyapunov exponent in position space is
$\gamma=\max\\{\gamma_{\varepsilon}^{+}(E),\gamma_{\varepsilon}^{-}(E)\\}$,
which is

$\begin{array}[]{r}\gamma(E)=\max\left\\{\operatorname{arcosh}\frac{|E+V+2|+|E+V-2|}{4},\right.\\\
\left.\operatorname{arcosh}\frac{|E-V+2|+|E-V-2|}{4}\right\\}.\end{array}$
(S6)

##  S-2 Derivation of Lyapunov exponent $\gamma_{m}$ in momentum space

In the main text, utilizing Fourier transform, the initial wave function
solution has been obtained. Relying on Jensen’s formula Jensen , our
calculation supposes that $f$ is an analytic function,
$a_{1},a_{2},\ldots,a_{n}$ are the zeros of $f$ in the interior of the unit
disc of the complex plane, and $f(0)\neq 0$. Then, we have the following
equality

$\ln|f|=\sum_{k=1}^{n}\ln\left({\left|a_{k}\right|}\right)+\frac{1}{2\pi}\int_{0}^{2\pi}\ln\left|f\left(e^{i\theta}\right)\right|d\theta.$
(S7)

Combining the initial wave function solution, the expression of Lyapunov
exponent can be written as

$\gamma_{m}(E)=\lim_{k\rightarrow\infty}\frac{1}{k-k_{0}}\ln\left|\frac{f_{k}}{f_{k_{0}}}\right|=\int\ln
g^{(1)}-\ln g^{(2)}d\theta,$

where $g^{(1)}=|-2\cos(2\pi\theta)+V+E|$, $g^{(2)}=\mid 2\cos(2\pi\theta)+$
$V-E\mid$.

Then, the first term of the rightmost side of $\gamma_{m}$ can be written as

$\displaystyle\int_{0}^{1}\ln g^{(1)}d\theta$
$\displaystyle=\int_{0}^{1}\ln|-2\cos(2\pi\theta)+V+E|d\theta$ (S8)
$\displaystyle=\int_{0}^{1}\ln|-e^{i2\pi\theta}-e^{-i2\pi\theta}+V+E|d\theta$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\ln|e^{i2\theta}+1-(V+E)e^{i\theta}|d\theta.$

Applying Jensen’s formula Jensen , the integral calculation of Eq. S8 can be
transformed to the calculation of roots of Eq. S9 in the unit disc, let
$x=e^{i\theta}$,

$x^{2}+1-(V+E)x=0.$ (S9)

After some mathematical calculations, we obtain

$\left\\{\begin{array}[]{l}\text{ when }|V+E|<2,\int_{0}^{1}\ln
g^{(1)}d\theta=0,\\\ \text{ when }|V+E|>2,\int_{0}^{1}\ln
g^{(1)}d\theta=\ln\frac{|E+V|+\sqrt{(E+V)^{2}-4}}{2}.\end{array}\right.$ (S10)

Under the similar process, with regard to the second term of $\gamma_{m}$, we
also obtain

$\left\\{\begin{array}[]{l}\text{ when }|V-E|<2,\int_{0}^{1}\ln
g^{(2)}d\theta=0,\\\ \text{ when }|V-E|>2,\int_{0}^{1}\ln
g^{(2)}d\theta=\ln\frac{|E-V|+\sqrt{(E-V)^{2}-4}}{2}.\end{array}\right.$ (S11)

Interesting, when both $|V+E|<2$ and $|V-E|<2$ hold, we have

$\gamma_{m}(E)=\int\ln g^{(1)}-\ln g^{(2)}d\theta=0.$ (S12)

When both $|V+E|<2$ and $|V-E|>2$ hold, then

$\gamma_{m}(E)=\int\ln g^{(1)}-\ln
g^{(2)}d\theta=0-\ln\frac{|E-V|+\sqrt{(E-V)^{2}-4}}{2}<0.$ (S13)

When both $|V+E|>2$ and $|V-E|<2$ hold, then

$\gamma_{m}(E)=\int\ln g^{(1)}-\ln
g^{(2)}d\theta=\ln\frac{|E+V|+\sqrt{(E+V)^{2}-4}}{2}-0>0.$ (S14)

When both $|V+E|>2$ and $|V-E|>2$ hold, then

$\gamma_{m}(E)=\int\ln g^{(1)}-\ln
g^{(2)}d\theta=\ln\frac{|E+V|+\sqrt{(E+V)^{2}-4}}{2}-\ln\frac{|E-V|+\sqrt{(E-V)^{2}-4}}{2}\neq
0.$ (S15)

To sum up the above calculations, only when both $|V+E|<2$ and $|V-E|<2$ hold,
namely the eigenvalues $E\in[V-2,2-V]$, $\gamma_{m}$ is equal to 0.

##  S-3 More numerical verification

In this section, we provide more numerical validation to strengthen the
credibility of our theoretical prediction. We have calculated the wave
functions of different energy levels at the same size, and the same energy
level at different sizes. The numerical results are as expected, and the
corresponding eigenstates dispaly self-similarity. As shown in Fig. S1 and
Fig. S2, it clearly illusrates that the scaled smaller peaks are very similar
to the largest peak. More interesting, in Fig. S2, three left small peaks
satisfy the scale invariance with a certain multiple $5.7/3.2\approx
3.2/1.8\approx 1.78$; whereas three right small peaks satisfy the scale
invariance with a certain multiple $13.9/5.7\approx 34/13.9\approx 2.44$. The
different scaling factors indicate that there is more than one fractal
structure in the wave function, which exactly corresponds to the multifractal
theory in the introduction.

Figure S1: (Color online) The black curve represents wave function of $E=0.5$
with the parameter $V=1$. The red curves represent three wave function peaks
after magnifying. It clearly shows that the scaled two smaller peaks are very
similar to the largest peak after twice magnification. The total number of
sites is set to be $L=6765$.
Figure S2: (Color online) The black curve represents wave function of $E=-1$
with the parameter $V=1$. The red curves represent three wave function peaks
after magnifying. The scaled multiples of three left small peaks are 5.7, 3.2
and 1.8 in turn; the scaled multiples of three right small peaks are 5.7, 13.9
and 34 in turn. The total number of sites is set to be $L=10946$.

## References

  * (1) A. Avila, Global theory of one-frequency Schrödinger operators, Acta. Math. 1, 215, (2015).
  * (2) J. Jensen, Sur un nouvel et important théorème de la théorie des fonctions, Acta Mathematica 1,22, (1899).


11institutetext: National Radio Astronomy Observatory, 520 Edgemont Rd.,
Charlottesville, VA 22903, USA 11email<EMAIL_ADDRESS>22institutetext:
Chester F. Carlson Center for Imaging Science, School of Physics & Astronomy,
and Laboratory for Multiwavelength Astrophysics, Rochester Institute of
Technology, 54 Lomb Memorial Drive, Rochester NY 14623 USA 33institutetext:
Institut Universitaire de France 44institutetext: Université Grenoble Alpes,
CNRS, IPAG, F-38000 Grenoble, France

# Tracing Molecular Stratification within an Edge-on Protoplanetary Disk

D. Ruíz-Rodríguez 1122 J. Kastner 22 P. Hily-Blant 33 T. Forveille 44

High-resolution observations of edge-on proto-planetary disks in emission from
molecular species sampling different critical densities and formation pathways
offer the opportunity to trace the vertical chemical and physical structures
of protoplanetary disks. Among the problems that can be addressed is the
origin and significance of the bright CN emission that is a ubiquitous feature
of disks. Based on analysis of sub-arcsecond resolution Atacama Large
Millimeter Array (ALMA) archival data for the edge-on Flying Saucer disk
(2MASS J16281370-2431391), we establish the vertical and radial
differentiation of the disk CN emitting regions with respect to those of 12CO
and CS, and we model the disk physical conditions from which the CN emission
arises. We demonstrate that the disk 12CO (2-1), CN (2-1), and CS J=5-4
emitting regions decrease in scale height above the midplane, such that 12CO,
CN, and CS trace layers of increasing density and decreasing temperature. We
find that at radii $>$ 100 au from the central star, CN emission arises
predominantly from intermediate layers, while in the inner region of the disk,
CN appears to arise from layers closer to the midplane. We investigate disk
physical conditions within the CN emitting regions, as well as the ranges of
CN excitation temperature and column density, via RADEX non-LTE modeling of
the three brightest CN hyperfine lines. Near the disk midplane, where we
derive densities $\rm n_{H_{2}}$ $\sim$107 cm-3 at relatively low $\rm
T_{kin}$ ($\sim$ 12 K), we find that CN is thermalized, while sub-thermal,
non-LTE conditions appear to obtain for CN emission from higher (intermediate)
disk layers. We consider whether and how the particular spatial location and
excitation conditions of CN emission from the Flying Saucer can be related to
CN production that is governed, radially and vertically, by the degree of
irradiation of the flared disk by X-rays and UV photons from the central star.

###### Key Words.:

Astrochemistry; Stars: pre-main sequence; Submillimeter: planetary systems.

##  1 Introduction

Detailed chemical studies of dusty circumstellar (protoplanetary) disks
orbiting young stars represent necessary steps toward understanding how
exoplanets and exocomets obtain their chemical inventories. The unprecedented
combination of sensitivity and spatial and spectral resolution afforded by the
Atacama Large Millimeter Array (ALMA) is now enabling such investigations
(e.g., Qi et al., 2013; Öberg et al., 2015; Guzmán et al., 2017; Hily-Blant et
al., 2017; Bergner et al., 2018; Kastner et al., 2018; Semenov et al., 2018;
Qi et al., 2019; Öberg & Wordsworth, 2019).

Among the problems that can be addressed via spatially resolved ALMA studies
of protoplanetary disk chemistry is the origin and significance of the bright
CN emission that is a ubiquitous feature of disks (e.g., Kastner et al., 1997,
2008, 2014; Thi et al., 2004; Sacco et al., 2014; Guilloteau et al., 2014).
Different astrochemical models of protoplanetary disk CN production, most of
which attempt to account for the influence of stellar high-energy irradiation,
have made a variety of predictions as to CN abundances and locations (e.g.
Lepp & Dalgarno, 1996; Stäuber et al., 2005; Chapillon et al., 2012; Cleeves
et al., 2013; Cazzoletti et al., 2018). In particular, most models predict
that CN is mostly produced within disk surface layers via UV irradiation
(either through photodissociation of HCN or via vibrational excitation of H2;
Cazzoletti et al., 2018; Visser et al., 2018). However, the low CN kinetic
temperatures inferred from observations ($T\leq 25$ K; e.g. Kastner et al.,
2014; Teague et al., 2016; Hily-Blant et al., 2017) suggest that disk CN
production may be driven by deeper-penetrating X-ray radiation from the
central T Tauri stars (Kastner et al., 2014).

Confirming or ruling out whether CN can be formed deep within protoplanetary
disks would have several interesting implications. Foremost, it would
contribute to establishing whether CN is a photodissociation product of HCN,
or is formed via other channels, such as reactions of N with C2H (Cazzoletti
et al., 2018). This, in turn, would help to constrain the physical and
chemical processes that took place during the formation of the various
isotopic reservoirs of nitrogen present in the Solar System (Guzmán et al.,
2017; Hily-Blant et al., 2017, 2019).

Despite the status of CN as a signature chemical tracer within disks — and the
foregoing open questions concerning its production and its utility as a tracer
of UV irradiation and nitrogen isotopic ratios — the only ALMA CN images of
disks published thus far are those for the nearly pole-on TW Hya and a handful
of moderately inclined disks (Teague et al., 2016; Hily-Blant et al., 2017;
van Terwisga et al., 2018; Teague & Loomis, 2020). The ring-like CN emission
morphologies revealed by these images motivated recent model calculations by
Cazzoletti et al. (2018) and, indeed, the observed morphologies are well
reproduced by these simulations. However, the degeneracy inherent to low-
inclination viewing geometries has left constraints on the scale height of CN
emission entirely in the realm of the aforementioned temperature-based
inferences.

To more directly investigate the vertical distribution of CN molecules within
protoplanetary disks, we have exploited archival ALMA data obtained for the
nearly edge-on Flying Saucer disk (2MASS J16281370$-$2431391). The CN data
analyzed here were obtained as part of the comprehensive study of the vertical
and radial physical (density and temperature) structure of that disk carried
out by Guilloteau et al. (2016), and subsequently by Dutrey et al. (2017).

The Flying Saucer system consists of a 0.57 M⊙ young stellar object surrounded
by a nearly edge-on disk ($i\approx 87^{\circ}$) of mass 3.5 $\times$ 10-5 M⊙.
Its distance is estimated at $\sim$ 120 pc (Loinard et al., 2008). As no
parallax is available for the Flying Saucer system in Gaia DR2 (Gaia
Collaboration et al., 2018), we adopt this distance for all calculations in
this paper. The aforementioned ALMA observations revealed the 12CO J = 2-1
emission lying well above and below the disk plane, arising from upper and
lower layers, extending to disk radii of at least 330 au. By contrast, the CS
J = 5-4 molecular emission is confined towards the midplane, and extends to
$\sim$ 300 au in radius (Dutrey et al., 2017). Contrary to the radial
extension of the CS and 12CO (hereafter CO) emission, the dust continuum is
limited to a more compact region with a radius of 190 au and, similarly to the
CS emission, is confined towards the midplane, where the dust temperature is
remarkably low ($T_{Dust}\sim$ 5-8 K at 100 au; Guilloteau et al., 2016).
Modeling of continuum and near-IR emission suggests that the disk is in an
advanced state of dust settling, with scale heights of 22.5 $\pm$ 1.5 au for
small grains ($\sim$ 1 $\mu$m) and 12.7 $\pm$ 0.3 au for large grains ($\sim$
1 mm) at 100 au (Grosso et al., 2003).

Indeed, CO, CS, and 242 GHz continuum emission reveal a vertical
stratification of the dust and gas in the system (Dutrey et al., 2017). These
stratified disk density and temperature tracers hence potentially afford a
direct means to constrain CN production pathways under different thermal and
density conditions. To exploit this potential, we have analyzed the CN
emission line data for the Flying Saucer disk that were obtained during the
same ALMA observations that yielded the maps of CO, CS, and continuum (Dutrey
et al., 2017). Sections 2 and 3 present details of the observations, and the
CS, CN, and CO molecular emission line maps resulting from our reprocessing of
the archival ALMA data. In Sec. 3, we present an analysis of the disk vertical
intensity distributions and CN emission line profiles, together with the
derivation of the disk physical parameters. Section 4 presents a discussion of
the disk structure, CN abundances, and potential triggers of CN production in
the disk at different scale heights and radii.

Table 1: Detected CN (N = 2-1) HF transitions.

HF Componenta | Transition | $\nu$ |  Aul |  gu |  Relative Intensitiesb
---|---|---|---|---|---
|  | [GHz] |  |  |
A |  J = 5/2 $\rightarrow$ 3/2, F = 5/2 $\rightarrow$ 3/2 | 226.8741908 |  9.62$\times$10-5 | 6 | 0.1805
B |  J = 5/2 $\rightarrow$ 3/2, F = 7/2 $\rightarrow$ 5/2 | 226.8747813 |  1.14$\times$10-4 | 8 | 0.2860
C |  J = 5/2 $\rightarrow$ 3/2, F = 3/2 $\rightarrow$ 1/2 | 226.8758960 |  8.59$\times$10-5 | 4 | 0.1074
D |  J = 5/2 $\rightarrow$ 3/2, F = 3/2 $\rightarrow$ 3/2 | 226.8874202 |  2.73$\times$10-5 | 4 | 0.0342
E |  J = 5/2 $\rightarrow$ 3/2, F = 5/2 $\rightarrow$ 5/2 | 226.8921280 |  1.81$\times$10-5 | 6 | 0.0340

a CN (2-1) has 19 HF components. We list here the 5 brightest (hence well-
detected) HF components within the observed bandpass, referred to throughout
this paper as A-E. Molecular data from Skatrud et al. (1983).

b Theoretical relative intensities under LTE conditions of the five brightest
HF transitions of CN within the bandwidth, adapted from Table 2 of Punzi et
al. (2015).

##  2 Observations

Figure 1: (a) Integrated emission map of CS from -0.9 to 8.6 km s-1 with
contours showing 10, 20, 30, 40, 50, and 60 $\sigma$ emission ($\sigma$ = 55
mJy beam-1 km s-1); the 0.53${}^{{}^{\prime\prime}}$$\times$
0.51${}^{{}^{\prime\prime}}$ at P.A. -57.8o beam is indicated in the bottom
left corner. (b) Integrated emission map of CN from 0.8 to 10.4 km s-1 with
contours showing 10, 20, 30, 40, 50, 60, 70, and 80 $\sigma$ emission $\sigma$
= 60 mJy beam-1 km s-1, and a 0.57${}^{{}^{\prime\prime}}$$\times$
0.55${}^{{}^{\prime\prime}}$ at P.A. -55.0o beam. (c) Integrated emission map
of CO from -0.9 to 8.6 km s-1 with contours showing 10, 20, 30, and 40
$\sigma$ emission ($\sigma$ = 90 mJy beam-1 km s-1), and a
0.56${}^{{}^{\prime\prime}}$$\times$ 0.54${}^{{}^{\prime\prime}}$ at P.A.
-63.4o beam. (d) ALMA 1.236 mm dust continuum image (natural weighting) with a
0.64${}^{{}^{\prime\prime}}$$\times$ 0.60${}^{{}^{\prime\prime}}$ at P.A. -56o
beam. At the bottom of each panel are shown the radial intensity profiles of
the integrated CS, CN and CO emissions and 242 GHz continuum. The radial
profiles are extracted along a cut through the disk midplane (P.A.= 3o) and
displayed by cyan lines in each image. Similarly, at the right of each panel
are shown the representative vertical intensity profiles of the integrated CS,
CN and CO emissions. The vertical profiles are extracted along a cut at a
radius of 170 au for CS, CN, and CO, and at 120 au for the 242 GHz continuum.
See Section 3.1.

The CN N=2-1, CS J = 5-4 and CO J = 2-1 observations of the Flying Saucer
presented in this work were obtained on May 17 and 18, 2015 as part of the
Cycle 2 program 2013.1.00387.S (PI: Guilloteau). Using 37 twelve meter
antennas, the range of unprojected baseline length is from 15 m to 560 m, and
the total observation time is $\sim$ 55 minutes. As previously noted, a
detailed description and analysis of the continuum, CO, and CS data was
presented in Guilloteau et al. (2016) and Dutrey et al. (2017).

Here, we reprocessed the ALMA data. To that end, the ALMA visibility data were
edited, calibrated and imaged using the pipeline version r39732 in CASA 4.2.2.
To image the data, we used the TCLEAN algorithm. We performed self-calibration
on the source continuum emission with three rounds of phase calibration. To
achieve a good balance between sensitivity and angular resolution the Briggs
weighting parameter $R$ was set to 0.5 for the continuum map and spectral line
image cubes. For the continuum map this yields a synthesized beam of
0.53${}^{{}^{\prime\prime}}$$\times$ 0.52${}^{{}^{\prime\prime}}$ at P.A.
-65.1o. The peak 242 GHz continuum flux is 8.5 mJy beam-1 and the r.m.s. is
0.2 mJy beam-1. Given the relatively weak continuum emission, continuum
subtraction should not remove significant line flux (Weaver et al., 2018). For
the lines of interest, we then subtracted the continuum emission by fitting a
first-order polynomial to the continuum in the $uv$-plane, and produced
spectral data cubes with 0.06 km s-1 velocity resolution. The image cubes were
constructed on a 256 $\times$ 256 pixel grid with 0.1′′ pixel size. For CS,
the resulting rms was $\sim$2 mJy beam-1 ($\sim$ 0.10 K) with a synthesized
beam of 0.53${}^{{}^{\prime\prime}}$$\times$ 0.51${}^{{}^{\prime\prime}}$
($\sim$ 60 au) at P.A. -57.8o. For CO, the rms is $\sim$2.5 mJy beam-1 ($\sim$
0.14 K) with a synthesized beam of 0.56${}^{{}^{\prime\prime}}$$\times$
0.54${}^{{}^{\prime\prime}}$ ($\sim$65 au) at P.A. -63.4o. For CN, the rms was
$\sim$2 mJy beam-1 ($\sim$ 0.11 K) with a synthesized beam of
0.57${}^{{}^{\prime\prime}}$$\times$ 0.55${}^{{}^{\prime\prime}}$ ($\sim$65
au) at P.A. -55.0o. In the latter case, the five brightest (of nineteen) CN
(2-1) hyperfine (HF) components are detected within the frequency range
covered by these observations, 226.87 - 226.89 GHz. These CN HF lines at
226.8741908, 226.8747813, 226.8758960, 226.8874202, 226.8921280 GHz are
labeled as ‘A’, ‘B’, ‘C’, ‘D’ and ‘E’ in Table 1, respectively.

We constructed moment 0 maps from each of the three spectral line data cubes
(Sec. 3.1), including a map merging CN components A–E. by clipping the signal
at $<$1$\sigma_{rms}$. The total continuum flux is 37.4 $\pm$ 0.88 mJy at 242
GHz, which is in good agreement with 35 mJy with $\sim$7$\%$ calibration
uncertainty obtained by Dutrey et al. (2017) from the same ALMA dataset. The
integrated emission line intensities over the entire emitting region are 1.60
$\pm$ 0.05 Jy km s-1, 3.90 $\pm$ 0.12 Jy km s-1, and 1.80 $\pm$ 0.16 Jy km s-1
for CS, CN, and CO, respectively. These flux uncertainties are estimated with
the sum in quadrature of $\lesssim$10$\%$ flux and rms (per the ALMA Cycle 2
Technical Handbook111https://almascience.eso.org/documents-and-
tools/cycle-2/alma-technical-handbook/). The CN line flux is consistent with
that previously obtained by Reboussin et al. (2015), 3.4 $\pm$ 0.5 Jy km s-1,
with the IRAM 30 m telescope.

##  3 Results

###  3.1 Integrated Intensity Maps and Scale Heights

Table 2: Results of vertical profile analysis.

Data | Radius | FWHMa | Beam | Hrc
---|---|---|---|---
| [au] | [au] | [au] | [au]
242 GHz |  | 99.4 $\pm$ 0.3 | 63.8 | 32.0 $\pm$ 1.2
CS |  | 129.5 $\pm$ 0.9 | 63.3 | 48.0 $\pm$ 0.7
CN | 60 | 154.9 $\pm$ 0.7 | 68.2 | 59.1 $\pm$ 0.6
CO |  | 193.2 $\pm$ 1.7 | 67.5 | 76.9 $\pm$ 0.9
242 GHz |  | 99.7 $\pm$ 0.3 | 63.8 | 32.3 $\pm$ 1.1
CS |  | 143.5 $\pm$ 1.8 | 63.3 | 54.7 $\pm$ 0.9
CN | 100 | 169.1 $\pm$ 1.5 | 68.2 | 65.7 $\pm$ 0.8
CO |  | 127.6 $\pm$ 3.0b | 67.5 | 107.2 $\pm$ 1.6
242 GHz |  | 103.4 $\pm$ 0.7 | 63.8 | 38.1 $\pm$ 1.2
CS |  | 159.0 $\pm$ 4.8 | 63.3 | 61.9 $\pm$ 2.3
CN | 180 | 205.7 $\pm$ 4.5 | 68.2 | 82.4 $\pm$ 2.0
CO |  | 123.6 $\pm$ 3.9b | 67.5 | 138.3$\pm$ 2.2

a Parameters obtained from Gaussian fit. For details, see Sec. 3.1.

b Multi-gaussian fitting performed; see text.

c Characteristic emission height obtained after deconvolving FWHM from
beamsize.

Figure 2: Vertical scale heights estimated at various radii from the CO (blue
hexagons), CN (black diamonds), CS (red dots), and 242 GHz continuum (cyan
squares) integrated emission images. $H_{r}$ values are only measured over the
area where the signal is $>3\sigma$, while error bars display the standard
error of the mean.

Figure 1 shows the integrated intensity (moment 0) maps of the Flying Saucer
disk in the CS J = 5-4, CN N=2-1, and CO J = 2-1 emission lines together with
the 242 GHz continuum emission. Figure 1 also displays CS, CN and CO radial
intensity profiles along the midplane of the nearly edge-on disk in the images
and representative vertical intensity profiles extracted along a cut at 170 au
from the central source. In the CN moment 0 image and radial intensity
profile, the emission displays a deficit toward the central star, suggestive
of an edge-on ring structure. The CS moment 0 map and its radial profile
display a very similar, double-peaked morphology (see also Dutrey et al.,
2017). The western regions of the CO moment 0 image are adversely affected by
contaminating background emission (Dutrey et al., 2017), such that the radial
profile of CO emission is difficult to interpret. Emission from both CN and CS
peaks at radii of $\sim$70 au.

From the radial profiles, we find that the bulk of the CS line emission
($>$75$\%$) originates within approximately $\sim$200 au from the central
source and has a total extension of $\sim$300 au, as also reported by Dutrey
et al. (2017). By contrast, the CN emitting region extends to at least 360 au,
well beyond the maximum radius of CS emission. Moreover, at distances greater
than 180 au, where the CO emission intensity in the mid-plane drops abruptly
on the (uncontaminated) eastern side of the disk, the CN emission in the
moment 0 map remains bright and confined closer to the mid-plane.

To describe the molecular emission structure of the disk, we estimate the
characteristic vertical disk height of CS, CN, and CO emission ($H_{r}$) at
various radii from 60 au to $\sim$300 au. These vertical profiles were
extracted along a position angle (PA) of $\sim$3.6∘ measured east from north
(Dutrey et al., 2017). For the continuum, CS, and CN, which are marginally
resolved we use the FWHMs of Gaussian fits to the observed vertical profiles
to quantify $H_{r}$.

Assuming that the characteristic emission height at a given radius $r$ is
similar on the western and eastern sides of the disk, the estimated FWHM
values are averaged. Figure 1a and 1b illustrate representative profiles along
vertical cuts at a radius of 170 au on each side, verifying their similarity.
To obtain $H_{r}$, we then deconvolve the beam width ($\theta_{FWHM}$) from
the resulting FWHM, i.e.,

$\rm H_{r}=\frac{\left(\sqrt{\rm FWHM^{2}-\theta_{\rm
FWHM}^{2}}\right)}{2\sqrt{\rm 2\>\textup{ln}2}}.$ (1)

Note that the resulting characteristic emission height $H_{r}$ is likely
directly related to the vertical scale height of the molecular species in
question. However, given the caveats, in particular, the emission from a given
species only directly traces the disk density structure and abundance
gradients if optically thin — it is important not to equate the results for
$H_{r}$ to disk scale height measurements.

In the specific case of CO emission at $r\geq$ 100, whose vertical extent is
clearly larger than CS and CO, the shapes of the vertical profiles change
considerably as a function of $r$, as a result of CO freeze-out processes and
self-absorption by colder gas along the disk midplane (Dutrey et al., 2017).
These effects are manifested in the form of a lack of CO emission from the
midplane, generating double-peaked vertical intensity profiles (see Figure
1c). Hence, the disk scale height cannot be characterized by the above
approximation (Eq. 1). Instead, the height of the disk surface can be
estimated by accounting for the resulting double-peaked vertical profiles (see
Figure 1c), whose peak positions ($\mu$) increase as disk radius increases in
a flared disk. Thus, the upper and lower disk layers, also labeled as $H_{r}$,
can be computed from

$\rm H_{r}=\frac{\left(\sqrt{\rm
FWHM^{2}-\theta_{FWHM}^{2}}\right)}{2\sqrt{2\>\textup{ln}2}}+\left|\mu\right|,$
(2)

with $\rm FWHM$ defined as

$\rm FWHM=2\sqrt{\rm 2\>\textup{ln}2}\sigma,$ (3)

where $\mu$, and $\sigma$ are the peak position and the width of the Gaussian
fitted, respectively. Here, $\mu$ accounts for an offset in the peak of the
distribution. We then adopt the average of these (upper and lower) $H_{r}$
values. Bearing in mind that CO is contaminated by background emission from
extended molecular clouds at velocities around $\sim$1.8 km s-1 and in the
range 6-7 km s-1 (Guilloteau et al., 2016), we only analyze vertical profiles
from the least contaminated region, i.e. the east side of the disk. Overall,
$H_{r}$ values are only measured over the area where the signal is $>3\sigma$.
Table 2 displays the resulting $H_{r}$ values at 60, 100 and 200 au for the
CS, CN and CO molecular lines. The tabulated uncertainties only reflect the
formal errors on the Gaussian fits to the intensity profiles 222We find very
similar values of $H_{r}$ for CS and CN if we adopt the double-Gaussian
fitting approach used for CO, i.e., if we use Equation 2 rather than Equation
1 to define $H_{r}$..

In Figure 2, we illustrate $H_{r}$ as a function of radius in the disk. It is
immediately apparent that the peak emission from the three molecular species
is stratified, with CO, CN, and CS lying increasingly close to the disk
midplane (as traced by the dust continuum emission). We discuss these results
in more detail in Sec 4.1.

###  3.2 CN Emission

Figure 3: Channel maps of the CN emission in the Flying Saucer from 0.60 to
6.40 km s-1 with the velocity labels referring to the B component. The
spectral resolution is 0.2 km s-1 and the rms per channel is 14 mJy beam-1.The
synthesized beam size is represented in the lower left panel. The frequency
range covers the ‘A’, ‘B, and ‘C’ HF components, which are J = 5/2
$\rightarrow$ 3/2, F = 5/2 $\rightarrow$ 3/2, J = 5/2 $\rightarrow$ 3/2, F =
7/2 $\rightarrow$ 5/2, J = 5/2 $\rightarrow$ 3/2, F = 3/2 $\rightarrow$ 1/2,
respectively, and listed in Table 1.

Figure 3 presents channel maps of CN emission, where the A, B, C HF emission
components, indicated by the circles, are spatially resolved and appear
according to their frequency-equivalent velocity. These CN channel maps were
generated with central velocities from $+$0.8 to $+$6.6 km s-1 with respect to
the brightest HF component (labeled B in Table 1).

As is the case for the CO emission from the Flying Saucer (Dutrey et al.,
2017), these CN channel maps reveal the vertical structure of a flared,
gaseous disk, with blue-shifted material to the west (W) and redshifted
material to the east (E). However, in these CN channel maps the HF splitting
of the CN molecule visually distorts the typical emission expected from
Keplerian rotation seen at high inclinations. This visual distortion is caused
by the overlapping of the HF components in velocity space at different
locations in the disk, as is apparent from the PV diagram presented in Fig. 4.
Furthermore, the disk shows nearly symmetric CN emission, with an emission
deficit in the midplane at radii between $\sim$ 80 and 200 au. The CN emission
from this region seems to trace upper and lower surfaces, which, in some iso-
velocity curves, reconnect at the outer region of the disk, forming an outer
wall. This emission feature –highlighted by the B component– is most evident
in channels at projected velocities of $\approx$ 2.4 and $\pm$ 0.2 km s-1 and
5.0 and $\pm$ 0.2 km s-1 ($\approx$ 1.5 km s-1 from the systemic velocity). A
similar pattern is detected in the weaker components, D and E, which fall
outside the range plotted in Fig. 3 (see Appendix A), and in the channel maps
of CS emission (Fig. 11).

####  3.2.1 Position-Velocity Diagram

Figure 4: Position-velocity diagram along the disk major axis of the CN
emission. Dashed curve is a Keplerian velocity profile for a stellar mass of
0.61 M⊙ at an inclination of 87o and a systemic velocity of 3.70 km s-1 LSRK.
The VLSRK scale is with respect to the brightest HF component (B) detected at
frequency 226.8747813 GHz. Labels A, B, and C indicate the brightest HF
transitions, see Table 1. Cyan ellipses correspond to the regions from which
we extracted the spectra analyzed in Sec. 3.2.2.

Figure 4 displays a PV diagram centered on the three brightest CN (2-1) HF
components (A, B, and C; Table 1), obtained from a cut along the major axis of
the disk (i.e. P.A $\sim$ 93.5 deg.). In these as well as the fainter HF
components (D and E; see Appendix A), the PV diagram reveals a disk consistent
with Keplerian rotation – with redshifted emission to the east and blueshited
to the west – as is the case for the CO and CS lines (Dutrey et al., 2017).
Taking the brightest HF line (B component) as the reference frame, the PV
diagram is well reproduced by a Keplerian rotation profile in a geometrically
thin disk with an inclination of 87 $\pm$ 1.0 deg., a central mass of 0.60
$\pm$ 0.02 M⊙ and a Kinematic Local Standard of Rest velocity (vLSRK) of 3.70
$\pm$ 0.5 km s-1, assuming a distance of 120 pc. These values are in good
agreement with the previous estimate of 0.58$\pm$ 0.01 M⊙ obtained from a
model of the disk CS emission by Dutrey et al. (2017).

####  3.2.2 CN Line Profiles and Parameters

Figure 5: The observed CN spectra (red), extracted from the outer region of
the disk ($>$ 210 au), with the Gaussian fits (green). Top and bottom panels
show the HF structure extracted from the western and eastern sides of the
disk, respectively. Cyan lines in the observed spectra mark the brightest HF
components and have been scaled to match the peak intensity in the B component
of the observed spectrum. The dashed blue line represents the systemic
velocity of the system ($\sim$3.7 Km s-1). These data can be converted to Jy
km s-1 by using the unit conversion factor of $\sim$55 K.

The CN HF transitions are well resolved at radii larger than 210 au, on both
the east and west sides of the disk. Figure 5 presents the spectra extracted
from ellipse-shaped apertures, centered on regions spatially represented by
the cyan ellipses drawn in Fig. 4, for which the HF splitting is well
resolved. We used Gaussian fits to these line profiles to determine line
widths and relative integrated intensities of the five HF components A–E. The
results, for components A to C, are shown in Fig. 5 for the west and east
sides, and results for all 5 components are summarized in Table 3. Under
optically thin local thermodynamic equilibrium (LTE) conditions, the ratio of
the intensity of the HF components relative to the brightest component (A/B,
C/B, D/B, and E/B) would range between 0.1 and 0.7; the theoretical intensity
ratios of the HF components are listed in Table 3. Under optically thick
conditions and for a uniform excitation temperature, the lines would saturate,
reaching equal intensities in the limit of very large optical depth.

At radii $>$ 210 au, the A, B, C, D and E components are at most marginally
optically thick, as determined by the measured HF transition ratios shown in
Table 3. We also note that the estimated HF transition ratios within 210 au
(not listed) trend toward 1.0, suggesting the emission is becoming optically
thicker. Overall, the D and E components in the CN spectrum are smeared out
towards the central source, becoming unresolvable, while the emission in A and
C are boosted to intensities nearly comparable to the B component until those
components (A+B+C) are blended together into a broad profile centered close to
the systemic velocity (3.70 km s-1).

Table 3: Flux ratios for integrated spectra at radii $>$ 210 au

|  |  | Westb |  |  |  | Eastb |
---|---|---|---|---|---|---|---|---
HF Component |  Intensity Ratioa | Intensity Ratio | FHWM |  $\Delta$Vc |  | Intensity Ratio | FHWM |  $\rm\Delta$Vc
| (Theoretical) |  |  |  km s-1 |  |  |  km s-1 |
A | 0.6311 |  0.68 $\pm$ 0.05 |  0.40 $\pm$ 0.02 |  0.24 $\pm$ 0.01 |  |  0.54 $\pm$ 0.02 |  0.41 $\pm$ 0.01 |  0.25 $\pm$ 0.01
B | 1 | 1 |  0.41 $\pm$ 0.01 |  0.25 $\pm$ 0.00 |  | 1 |  0.51 $\pm$ 0.01 |  0.31 $\pm$ 0.01
C | 0.3755 |  0.47 $\pm$ 0.04 |  0.38 $\pm$ 0.02 |  0.23 $\pm$ 0.01 |  |  0.38 $\pm$ 0.02 |  0.47 $\pm$ 0.02 |  0.28 $\pm$ 0.01
D | 0.1196 |  0.15 $\pm$ 0.04 |  0.33 $\pm$ 0.06 |  0.20 $\pm$ 0.03 |  |  0.17 $\pm$ 0.02 |  0.42 $\pm$ 0.03 |  0.25 $\pm$ 0.02
E | 0.1189 |  0.10 $\pm$ 0.03 |  0.30 $\pm$ 0.08 |  0.18 $\pm$ 0.05 |  |  0.16 $\pm$ 0.02 |  0.40 $\pm$ 0.03 |  0.24 $\pm$ 0.02

a Theoretical ratio of intensities of HF transitions to the brightest HF
transition (J = 5/2 $\rightarrow$ 3/2, F = 7/2 $\rightarrow$ 5/2). The
relative intensities of HF transitions of CN are adapted from Table 2 of
(Punzi et al., 2015).

b Values obtained from Gaussian fitting (see Sec. 3.2.2).

C $\rm\Delta$V = $\frac{\rm FWHM}{\sqrt{4\>\textup{ln}2}}$.

###  3.3 CN Column Densities

The analysis of Sec. 3.2.2 provides indications of CN line opacities, but does
not address the key physical conditions that give rise to the emission, in
particular, CN column densities333Note that the $\rm N_{CN}$ values derived
here represent line-of-sight CN column densities, as opposed to CN column
densities measured vertically through the disk. ($\rm N_{CN}$) and excitation
temperatures ($\rm T_{ex}$). To ascertain the range of these and other
parameters (e.g., optical depth, $\rm\tau$, and H2 number density, $\rm
n_{H_{2}}$) as functions of position within the disk, we have modeled the CN
HF line structure using a grid-based parameter exploration approach that
exploits the non-LTE radiative transfer code RADEX under the so-called Large
Velocity Gradient (LVG) framework (van der Tak et al., 2007). RADEX uses an
escape probability formalism, with five parameters as inputs: $\rm n_{H_{2}}$,
$\rm N_{CN}$, kinetic temperature ($\rm T_{kin}$), background temperature, and
linewidth. To explore a range of thermal parameters consistent with previous
studies (e.g. Guilloteau et al., 2016; Dutrey et al., 2017), we constructed a
3-D grid consisting of 21 values of H2 density in the range $\rm n_{H_{2}}$ =
104 \- 109 cm-3, 21 values of CN column density in the range $\rm N_{CN}$ =
1011 \- 1017 cm-2, and 51 values of $\rm T_{kin}$ from 5 to 60 K. As part of
the fitting process (see below), we upsampled these RADEX 3D-grid models by
interpolating between grid points to increase the resolution of the parameter
space explored. The FWHM line width of all three HF components was held fixed
at 0.4 km s-1, which is the average of the individual line-widths derived in
Sec. 3.2.2. For the background temperature, we adopt 2.73 K. For CN emission,
RADEX then predicts $\rm T_{ex}$ and $\rm\tau$ values for the ‘A’, ‘B’, and
‘C’ components from each set of these input parameters ($\rm T_{kin}$, $\rm
n_{H_{2}}$, linewidth, $\rm N_{CN}$) using collision rates from Hily-Blant et
al. (2013). Since RADEX does not treat line overlaps, we restricted the
spectral modeling to disk regions where these CN HF structures are well
resolved (see below).

To perform HF line fitting covering the full spectral range of the ‘A’ $+$ ‘B’
$+$ ‘C’ components, we use a forward-fitting approach built into the PySpecKit
Python package (Ginsburg & Mirocha, 2011). For each grid point (i.e.,
combination of $\rm T_{kin}$, $\rm n_{H_{2}}$, and $\rm N_{CN}$), the RADEX-
computed $\rm T_{ex}$ and $\rm\tau$ parameters are used to generate synthetic
line profiles for each of the three HF components from

$\rm T_{MB}(\nu)=\rm
T_{o}\left(\frac{1}{e^{\frac{T_{o}}{T_{ex}}}-1}-\frac{1}{e^{\frac{T_{o}}{2.73}}-1}\right)\left(1-e^{-\tau}\right),$
(4)

where $\rm T_{MB}(\nu)$ is the predicted brightness temperature, $\rm
T_{o}(\nu)=\frac{h\nu}{k}$, $k$ is Boltzmann’s constant, and $h$ is Planck’s
constant (e.g. Mangum & Shirley, 2015). The complete synthetic spectrum of the
CN 2–1 lines for a given parameter set is obtained as the sum of intensities
of the individual HF components, each described by Eq. 4. This procedure is
validated a posteriori because the excitation temperatures of each component
are very similar and line overlap for the A and B components is modest (see
Figures 5 and 7). The set of synthetic spectra is then compared to the
observed spectra using a $\chi^{2}$ minimization process implemented in the
PySpecKit Python package, which employs the Levenberg-Marquardt optimization
method via the MPFIT algorithm (Markwardt, 2009, and references within).

The preceding analysis was carried out for eight specific outer disk regions
(r $>$ 210 au) where the HF components are spectrally well resolved, on both
the east and west sides of the disk; see Figure 6 and Table 4. Figure 7
displays the CN spectra extracted from these eight regions — each of diameter
one beamwidth ($\sim$0.6${}^{{}^{\prime\prime}}$), centered at upper and lower
intermediate layers, midplane, and outer midplane, on both the east and west
sides of the disk (Figure 6) — overlaid with the best-fit LVG models, i.e.,
the model with the lowest $\chi^{2}$ value. The values of $\rm N_{CN}$, $\rm
n_{H_{2}}$, $\rm T_{kin}$, $\rm T_{ex}$, and $\rm\tau$ obtained for each
region from this $\chi^{2}$ minimization process are listed in Table 4. We
caution that — as is typical in such modeling — there is an inherent
degeneracy, wherein the HF lines are well reproduced either by very dense and
cold gas or by lower density, warmer gas.

These results (Table 4) indicate that the gas density traced by CN in the
regions we modeled (Figure 6) is a few $\times 10^{6}$ cm-3, and that the gas
is quite cold, with a best-fit $\rm T_{kin}$ in the range 10–12 K. We
consistently find $\rm T_{ex}$ $<$ $\rm T_{kin}$, where $\rm T_{ex}$/$\rm
T_{kin}$ decreases as $\rm T_{kin}$ increases, suggesting the gas may be
somewhat subthermally excited. For $\rm T_{kin}$ $>$ 25 K, we find the models
fail to reproduce (underestimate) the intensity of the C component of the CN
HF structure. The best-fit CN column densities are consistently $\sim$1013
cm-2. The modeling results also indicate that the west side of the disk may be
somewhat colder. We caution, however, that the results for $\rm T_{kin}$ are
subject to systematic uncertainties (arising from, e.g., the various model
assumptions and the specific approach to $\chi^{2}$ minimization) that are
difficult to quantify. The possibility of such an east-west asymmetry in the
disk warrants followup via CN mapping at higher resolution.

Figure 6: Integrated intensity map of the CN emission overlaid with eight
enumerated circles representing the eight single-beamwidth
(0.3${}^{{}^{\prime\prime}}$ radius) regions from which we extract spectra at
each side of the disk at radii $>$ 210 au (see Table 4, Section 3.3, and Fig.
7). The synthesized beam size is indicated by the white ellipse. Vertical
green dashed lines are placed at radii of $\sim$210 au from the central star;
the horizontal green dashed line indicates the disk midplane..

Figure 7: Spectra of the integrated CN (2-1) emission from the upper and lower
surfaces, midplane and outer midplane of the disk along with the simultaneous
fits to components ‘A’, ‘B’ , and ‘C’; see Table 3. The best-fit LVG models
from the $\chi^{2}$ minimization are displayed as a green line. The fitted
parameters are the total optical depth ($\rm\tau$; ‘A’ + ‘B’ + ‘C’),
excitation temperature ($\rm T_{ex}$), and on logarithmic scale, CN column
density ($\rm N_{CN}$) and particle density ($\rm n_{H_{2}}$); see Table 4.
Spectra in the left and right columns are extracted from the west and east
sides of the disk, respectively. The fit residuals are shown as a function of
frequency at the bottom of each panel. The shaded gray area shows the error
computed from a signal-free part of the spectra.  Table 4: RADEX Model
Parameters.

No. | Region | Location | $\rm T_{kin}$ | $\rm n_{H_{2}}$ | $\rm N_{CN}$ | $\rm T_{ex}$ |  $\rm\tau$a
---|---|---|---|---|---|---|---
|  |  $\Delta$ R.A. |  $\Delta$ Dec. | [K] |  $\times$ 106 [cm-3] |  $\times$ 1013 [cm-2] | [K] |
|  | [au] | [au] |  |  |  |  |
1. | W. Upper Intermediate Layer |  $+$ 230 |  $+$ 50 | 9 | 6.8 | 3.5 | 7.3 | 2.8
2. | W. Midplane |  $+$ 230 |  $+$ 0.0 | 10 | 50.0 | 6.0 | 9.5 | 3.5
3. | W. Lower Intermediate Layer |  $+$ 230 |  $-$ 50 | 9 | 8.0 | 4.4 | 7.7 | 3.5
4. | W. Outer Midplane |  $+$ 280 |  $+$ 0.0 | 7 | 7.3 | 4.2 | 5.8 | 3.8
5. | E. Upper Intermediate Layer |  $-$ 230 |  $+$ 50 | 11 | 5.4 | 3.9 | 8.1 | 3.0
6. | E. Midplane |  $-$ 230 |  $+$ 0.0 | 12 | 74.0 | 5.7 | 11.6 | 2.8
7. | E. Lower Intermediate Layer |  $-$ 230 |  $-$ 50 | 11 | 3.2 | 6.4 | 7.4 | 4.9
8. | E. Outer Midplane |  $-$ 280 |  $+$ 0.0 | 8 | 7.0 | 3.7 | 6.5 | 3.1

a Sum of the opacities at line center of A, B and C HF components.

##  4 Discussion

Figure 8: Top: Three-color (RGB) image generated from CS, CN, and CO peak
intensity (CASA moment 8) maps. CS, CN, and CO are assigned to the red, green,
and blue channels, respectively; each color channel is scaled linearly.
Bottom: Schematic illustration of the different layers traced by CS
(midplane), CN (midplane to intermediate layers) and CO (upper disk surface)
emission.

###  4.1 Vertical Molecular Emission Structure

The analysis presented in Sec. 3.1 (see, in particular Fig. 2), establishes
that the CS, CN, and CO emission regions in the Flying Saucer disk are
vertically stratified. To further demonstrate this vertical stratification, we
have generated maps of the peak intensities of CN, CS, and CO (i.e., CASA
“moment 8” maps) from their respective image cubes (presented in Figs. 3, 11,
and 12). These peak emission images are rendered as an RGB image in Fig. 8a,
where they are compared with a schematic illustration of the disk structure in
these three tracers (Fig. 8b).

The structure illustrated in Fig. 8 can be summarized as follows: CO reaches
the largest vertical extension, delineating the disk surface, with a deficit
of CO near the midplane; CS emission is confined to the layers closest to the
cold disk midplane; and CN traces layers intermediate between these extremes,
albeit with some overlap with CS, in that the CN also appears to extend to
near the midplane. As discussed in Dutrey et al. (2017), the midplane deficit
of CO seems to correspond to the vertical settling of dust grains toward the
midplane as traced by the 242 GHz continuum emission (Fig. 1d), which forms an
elongated structure of about $\sim$190 au in radius (see also Guilloteau et
al., 2016). Of the three molecular tracers in the ALMA data thus far available
for the Flying Saucer, it is evident from Fig. 2 and Fig. 8a that CS most
closely traces this (vertically and radially confined) continuum emission,
although CS is clearly more extended (in both dimensions). The CN originates
from a layer that is radially more extended than and elevated with respect to
CS. In addition, there appears to be a dearth of CN emission at the midplane
at radii between $\sim$80 and 200 au; this is apparent in the CN velocity
channel maps (Fig. 3). At the outer regions of the disk, the CN emission
arises from the upper and lower intermediate layers and extends down to the
midplane, outlining the disk radial extension. Possible mechanisms that might
enhance CN production in these (intermediate) disk layers are discussed in
Sec. 4.3.

Teague et al. (2016) previously suggested such a vertically stratified
emission structure for CS, CN, and CO emission lines in the case of the nearly
face-on TW Hydrae disk (TW Hya). From analysis of multiple transitions of the
three species, these investigators derived CO, CN, and CS excitation
temperatures of $\sim$35, 25, and 12 K, respectively (see also Teague &
Loomis, 2020). On this basis, Teague et al. (2016) deduced that the emission
from the three molecules arises from decreasing heights above the disk
midplane. This conclusion is consistent with the data presently available for
the nearly edge-on Flying Saucer (Fig. 8), although we deduce colder
excitation (and kinetic) temperatures for CN in this disk (see Sec. 4.3).
However, the implied discrepancy in the thermal structures of the TW Hya and
Flying Saucer disks is perhaps not surprising, given the different stellar
parameters and ages of both systems. In particular, the difference in the
derived values of $\rm T_{kin}$ might lie with the fact that the TW Hya is a
late-K star of rather advanced age ($\sim$8 Myr; e.g. Donaldson et al., 2016),
and is a luminous source of high-energy photons, with LX = 1.4 $\times$ 1030
ergs s-1 (Kastner et al., 2002), whereas the Flying Saucer central star is a
younger, M-type protostar that has thus far gone undetected in X-rays
444Chandra and XMM-Newton X-ray observations of a $\rho$ Oph field that
included Flying Saucer were obtained in 2000 and 2003 with exposure times of
$\sim$ 5 ks and $\sim$ 34 ks, respectively. Neither observation yield a
detection at the position of the object. Given the disk’s edge-on orientation,
however, it is difficult to draw conclusions as to the intrinsic X-ray
luminosity of the central star.. Furthermore, the TW Hya disk has a CN
emission extension of $\sim$180 au (Vlemmings et al., 2019; Teague & Loomis,
2020), a factor of $\sim$2 smaller than that of the Flying Saucer disk.

###  4.2 CN Physical Conditions

The results of our RADEX analysis indicate that the CN traces cold gas, and is
optically thick, in the disk regions modeled via RADEX, i.e., $\rm T_{kin}$
$\sim$ 7-10 K (with total $\rm\tau$ $\sim$ 3-4) and $\rm T_{kin}$ $\sim$ 8-12
K (with total $\rm\tau$ $\sim$ 3-5) at radii $>$210 au and heights $<$60 au in
the western and eastern sides of the disk, respectively (Sec 3.3). The low
midplane temperatures appear compatible with a model of the temperature
structure of the Flying Saucer disk obtained by Dutrey et al. (2017) from
numerical simulations of the disk’s CO emission. Specifically, those authors
found $T\sim 5-10$ K in the disk midplane at offsets of $\sim$230-280 au
(Dutrey et al., 2017, their Fig. 8g), corresponding to the western and eastern
midplane regions we included in our analysis (Table 4). Furthermore, the
densities we obtain from our line profile modeling ($\sim$$10^{6}$ to $10^{7}$
cm-3; Table 4) are generally consistent with, if somewhat smaller than, the
densities predicted for these same (external) disk regions in the Dutrey et
al. (2017) modeling (their Fig. 7). Their modeling also indicates, however,
that $T_{kin}\sim$ 15–25 K at heights of $\sim$50 au at the same
($\sim$230-280 au) radial offsets. These disk regions correspond to the
“intermediate layers” for which we find $\rm T_{kin}$ $\sim$ 9–11 K based on
RADEX line profile fitting (Table 4), and are not significantly different from
the corresponding values in the midplane. This lack of a significant vertical
gradient in $\rm T_{kin}$ at radii $>$210 au, as inferred from CN emission,
may indicate that the standard protoplanetary disk model, on which the Dutrey
et al. (2017) modeling is based, may not well describe the vertical
temperature structure in the Flying Saucer.

On the other hand, our RADEX modeling also yields evidence that the CN
molecules are subthermally excited ($T_{ex}<T_{kin}$) away from the midplane,
with the divergence between $\rm T_{ex}$ and $\rm T_{kin}$ more evident on the
eastern side of the disk (Table 4). We caution that, as noted earlier, it
remains to be determined whether the apparent differences in the CN
thermalization inferred for the two sides of the disk reflect a real asymmetry
in disk physical conditions. Indeed, one difficulty with the low $\rm T_{kin}$
inferred here is that one would expect CN and CS to be severely depleted from
the gas phase in disk regions in which $\rm T_{kin}$ drops well below the
freeze-out temperature of CO (i.e., $\sim$17 - 19 K). However, subthermal
excitation becomes more likely in disk regions with a significant degree of
depletion of dust grains, as may characterize the Flying Saucer at radii
$>210$ au. Previous studies have shown that when the gas-to-dust ratio drops
to $\sim$20 (or $n_{H_{2}}$ drops by a factor $\sim$5), CN would be produced
under subthermal excitation conditions, although this is only expected in
regions where the density is of the order of a few 105 cm-3 (Chapillon et al.,
2012).

###  4.3 CN: Location and Formation

Cazzoletti et al. (2018) modeled CN emission from protoplanetary disks under
the assumption that CN formation mainly proceeds through far ultraviolet (FUV)
pumping of H2 to vibrationally excited levels, in addition to
photodissociation of HCN. As a main result, the highest CN abundance should be
located in the upper region of the disk, while in outer disk regions, the
emission extends down to the midplane, forming a CN emission ring. Although
the HF transitions are not considered specifically in this model, the CN
features observed in the Flying Saucer disk appear reasonably consistent with
the predicted formation pathway, where CN is produced near the surface and
outer regions of the disk (Fig. 3). However, it is apparent from Fig. 2 that
the CN layer lies between the CO and CS layers, clearly showing that CN
emission traces intermediate layers rather than disk surface layers. In
addition, the strong link between UV radiation and CN formation would appear
to require that most of the CN emission should arise from upper disk layers,
given that the FUV radiation field is attenuated by dust at moderate disk
depths (Gorti et al., 2009). Indeed, the notion that UV controls the abundance
of CN is consistent with the deficit of midplane CN emission at radii between
$\sim$80 and 200 au (see Fig. 3), matching the spatial extension of the 242
GHz continuum.

However, as noted earlier, we also find that the CN arises from very cold
($\rm T_{kin}$ $\lesssim$15 K) regions of the disk (Sec. 4.2). This suggests
that FUV radiation alone may be insufficient to explain the observed CN
abundances, and one must invoke more deeply penetrating sources of high energy
radiation, such as X-rays (see discussions in Kastner et al., 2014; Hily-Blant
et al., 2017). Indeed, the bulk of the CN intensity is localized within the
central $\sim$0.8${}^{{}^{\prime\prime}}$ ($\sim$ 100 au) in radius (Fig. 1b)
of the disk, a region where enhancement of CN production can be produced by
exposing a disk that is in advanced stages of settling to X-rays or a
combination of X-rays and FUV radiation (Stäuber et al., 2005, 2007).
Moreover, recent simulations from detailed physical/chemical disk models
indicate that the distribution of radicals such as CN could be present as
multiple-layered structures (e.g. Ruaud & Gorti, 2019). In fact, several
channel maps of the CN line, Fig. 3, appear to show emission arising from the
midplane, especially at the inner part of the CN emission ‘ring’, while there
are velocities in the channel maps at which it appears that CN arises from
intermediate layers, at larger disk radii. These structures are qualitatively
similar to the predictions of the Ruaud & Gorti (2019) models. Both emission
regions may contribute to low-temperature CN emission, such as we infer from
our RADEX modeling (Sec. 3.3). Higher-resolution ALMA imaging spectroscopy of
multiple transitions of CN would help establish whether there are indeed two
distinct disk layers involved in efficient CN production.

##  5 Summary

Using archival ALMA data, we have investigated the vertical chemical and
physical structure of the edge-on Flying Saucer disk. This is in order to
constrain CN formation pathways that took place during disk evolution and
planetary formation. In particular, we used RADEX models to generate CN (2-1)
synthetic spectra, and thus predict the disk thermal conditions that
characterize the observed CN emission from the Flying Saucer disk. As a
comparative approach, we have used molecular lines tracing different disk
physical conditions, specifically12CO and CS, in order to further constrain
spatially the location of the CN emission. Our main results are as follows:

  * •

We demonstrate that 12CO, CN, and CS molecular lines can be used to trace the
vertical structure of edge-on disks. We estimated the scale heights of these
molecules, and they reveal a stratified structure within the flared Flying
Saucer disk. Basically, 12CO traces the low-density disk surface, while CN and
CS trace the denser intermediate layers and midplane region, respectively.

  * •

The analysis of the CN channel maps reveals complex structures superimposed on
Keplerian gas kinematics of the Flying Saucer disk. An outstanding feature is
the CN emission deficit near the midplane at radii $>$ 100 au, predominantly
arising from intermediate layers. However, in the inner regions of disk, CN
emission appears to arise from layers closer to the midplane. This variation
of the CN vertical location along the disk might be related to the formation
route of CN, whose production is likely related to the degree of irradiation
of the flared disk by X-rays and UV photons from the central star.

  * •

Based on RADEX non-LTE modeling of the three brightest HF components of the
CN(2-1) line, we estimated thermal parameters in different regions of the
disk. Our results characterized a disk with thermal emission from the midplane
and sub-thermal conditions in intermediate layers. Near the disk midplane,
where we derive densities $\rm n_{H_{2}}$ $\sim$107 cm-3 at relatively low
$\rm T_{kin}$ ($\sim$ 12 K), we find that CN is thermalized, while sub-
thermal, non-LTE conditions appear to obtain for CN emission from higher
(intermediate) disk layers.

  * •

Higher-resolution ALMA imaging spectroscopy of multiple transitions of CN is
required to confirm the cold outer disk temperatures we infer from our
radiative transfer modeling. Such followup observations would also help
establish whether there exist multiple, distinct disk layers involved in
efficient CN production.

###### Acknowledgements.

We wish to thank the anonymous referee for providing insightful comments. This
paper makes use of the following ALMA data: ADS/JAO.ALMA$\\#$2013.1.00387.S.
ALMA is a partnership of ESO (representing its member states), NSF (USA) and
NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio
Astronomy Observatory is a facility of the National Science Foundation
operated under cooperative agreement by Associated Universities, Inc.

## References

  * Bergner et al. (2018) Bergner, J. B., Guzmán, V. G., Öberg, K. I., Loomis, R. A., & Pegues, J. 2018, ApJ, 857, 69
  * Cazzoletti et al. (2018) Cazzoletti, P., van Dishoeck, E. F., Visser, R., Facchini, S., & Bruderer, S. 2018, A&A, 609, A93
  * Chapillon et al. (2012) Chapillon, E., Guilloteau, S., Dutrey, A., Piétu, V., & Guélin, M. 2012, A&A, 537, A60
  * Cleeves et al. (2013) Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5
  * Donaldson et al. (2016) Donaldson, J. K., Weinberger, A. J., Gagné, J., et al. 2016, ApJ, 833, 95
  * Dutrey et al. (2017) Dutrey, A., Guilloteau, S., Piétu, V., et al. 2017, A&A, 607, A130
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, ArXiv e-prints [arXiv:1804.09365]
  * Ginsburg & Mirocha (2011) Ginsburg, A. & Mirocha, J. 2011, PySpecKit: Python Spectroscopic Toolkit
  * Gorti et al. (2009) Gorti, U., Dullemond, C. P., & Hollenbach, D. 2009, ApJ, 705, 1237
  * Grosso et al. (2003) Grosso, N., Alves, J., Wood, K., et al. 2003, ApJ, 586, 296
  * Guilloteau et al. (2016) Guilloteau, S., Piétu, V., Chapillon, E., et al. 2016, A&A, 586, L1
  * Guilloteau et al. (2014) Guilloteau, S., Simon, M., Piétu, V., et al. 2014, A&A, 567, A117
  * Guzmán et al. (2017) Guzmán, V. V., Öberg, K. I., Huang, J., Loomis, R., & Qi, C. 2017, ApJ, 836, 30
  * Hily-Blant et al. (2017) Hily-Blant, P., Magalhaes, V., Kastner, J., et al. 2017, A&A, 603, L6
  * Hily-Blant et al. (2019) Hily-Blant, P., Magalhaes de Souza, V., Kastner, J., & Forveille, T. 2019, A&A, 632, L12
  * Hily-Blant et al. (2013) Hily-Blant, P., Pineau des Forêts, G., Faure, A., Le Gal, R., & Padovani, M. 2013, A&A, 557, A65
  * Kastner et al. (2014) Kastner, J. H., Hily-Blant, P., Rodriguez, D. R., Punzi, K., & Forveille, T. 2014, ApJ, 793, 55
  * Kastner et al. (2002) Kastner, J. H., Huenemoerder, D. P., Schulz, N. S., Canizares, C. R., & Weintraub, D. A. 2002, ApJ, 567, 434
  * Kastner et al. (2018) Kastner, J. H., Qi, C., Dickson-Vandervelde, D. A., et al. 2018, ApJ, 863, 106
  * Kastner et al. (2008) Kastner, J. H., Zuckerman, B., & Bessell, M. 2008, A&A, 491, 829
  * Kastner et al. (1997) Kastner, J. H., Zuckerman, B., Weintraub, D. A., & Forveille, T. 1997, Science, 277, 67
  * Lepp & Dalgarno (1996) Lepp, S. & Dalgarno, A. 1996, A&A, 306, L21
  * Loinard et al. (2008) Loinard, L., Torres, R. M., Mioduszewski, A. J., & Rodríguez, L. F. 2008, ApJ, 675, L29
  * Mangum & Shirley (2015) Mangum, J. G. & Shirley, Y. L. 2015, PASP, 127, 266
  * Markwardt (2009) Markwardt, C. B. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 411, Astronomical Data Analysis Software and Systems XVIII, ed. D. A. Bohlender, D. Durand, & P. Dowler, 251
  * Öberg et al. (2015) Öberg, K. I., Guzmán, V. V., Furuya, K., et al. 2015, Nature, 520, 198
  * Öberg & Wordsworth (2019) Öberg, K. I. & Wordsworth, R. 2019, AJ, 158, 194
  * Punzi et al. (2015) Punzi, K. M., Hily-Blant, P., Kastner, J. H., Sacco, G. G., & Forveille, T. 2015, ApJ, 805, 147
  * Qi et al. (2019) Qi, C., Öberg, K. I., Espaillat, C. C., et al. 2019, ApJ, 882, 160
  * Qi et al. (2013) Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013, Science, 341, 630
  * Reboussin et al. (2015) Reboussin, L., Guilloteau, S., Simon, M., et al. 2015, A&A, 578, A31
  * Ruaud & Gorti (2019) Ruaud, M. & Gorti, U. 2019, ApJ, 885, 146
  * Sacco et al. (2014) Sacco, G. G., Kastner, J. H., Forveille, T., et al. 2014, A&A, 561, A42
  * Semenov et al. (2018) Semenov, D., Favre, C., Fedele, D., et al. 2018, A&A, 617, A28
  * Skatrud et al. (1983) Skatrud, D. D., de Lucia, F. C., Blake, G. A., & Sastry, K. V. L. N. 1983, Journal of Molecular Spectroscopy, 99, 35
  * Stäuber et al. (2007) Stäuber, P., Benz, A. O., Jørgensen, J. K., et al. 2007, A&A, 466, 977
  * Stäuber et al. (2005) Stäuber, P., Doty, S. D., van Dishoeck, E. F., & Benz, A. O. 2005, A&A, 440, 949
  * Teague et al. (2016) Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49
  * Teague & Loomis (2020) Teague, R. & Loomis, R. 2020, ApJ, 899, 157
  * Thi et al. (2004) Thi, W., van Zadelhoff, G., & van Dishoeck, E. F. 2004, A&A, 425, 955
  * van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627
  * van Terwisga et al. (2018) van Terwisga, S. E., van Dishoeck, E. F., Cazzoletti, P., et al. 2018, arXiv e-prints [arXiv:1811.03071]
  * Visser et al. (2018) Visser, R., Bruderer, S., Cazzoletti, P., et al. 2018, A&A, 615, A75
  * Vlemmings et al. (2019) Vlemmings, W. H. T., Lankhaar, B., Cazzoletti, P., et al. 2019, A&A, 624, L7
  * Weaver et al. (2018) Weaver, E., Isella, A., & Boehler, Y. 2018, ApJ, 853, 113

##  Appendix A CN HF Components

Figure 9: Channel maps of the CN emission in the Flying Saucer from -4.60 to
-0.40 km -1, and from 1.40 to 5.60 km -1. The frequency range is presented for
the ‘D’ (Top), and ‘E’ (Bottom) HF components, which are J = 5/2 $\rightarrow$
3/2, F = 3/2 $\rightarrow$ 3/2, and J = 5/2 $\rightarrow$ 3/2, F = 5/2
$\rightarrow$ 5/2, respectively (see Table 3). The spectral resolution is 0.3
km s-1 and the rms per channel is 2.9 mJy beam-1. The synthesized beam size is
represented in the lower left panel.

The CN N=(2-1) HF components channel maps of the C and D components are
displayed in Fig. 9. These satellite components also present a typical
Keplerian rotation in a disk. In Figure 10, we present PV diagrams along the
disk major axis overlaid with curves representing Keplerian rotation in a
geometrically thin disk with an inclination of 87 $\pm$ 1.0 deg, a central
mass of 0.60 $\pm$ 0.02 M⊙ and vLSRK of 3.70 $\pm$ 0.5 km s-1, assuming a
distance of 120 pc.

Figure 10: Position-velocity diagram along the disk major axis of the CN
emission. Dashed curve is a Keplerian velocity profile for a stellar mass of
0.61 M⊙ at an inclination of 87o and a systemic velocity of 3.70 km s-1 LSRK.
Labels A, B, C, D and E indicate the 5 observed HF transitions, see Table 3.

##  Appendix B CS and CO emission: Channel maps

Channel maps of CS and CO toward Flying Saucer disk are displayed in Figs. 11
and 12. The velocity range is similar to channel maps for CN presented in Fig.
3 for a direct comparison.

Figure 11: Channel maps of the CS emission in the Flying Saucer from 0.60 to
6.40 km s-1. The spectral resolution is 0.2 km s-1 and the rms per channel is
3.0 mJy beam-1.The synthesized beam size is represented in the lower left
panel.  Figure 12: Channel maps of the CO emission in the Flying Saucer from
0.60 to 6.40 km s-1. The spectral resolution is 0.2 km s-1 and the rms per
channel is 5.0 mJy beam-1.The synthesized beam size is represented in the
lower left panel.