\section{Introduction}
Let $(X,L)$ be a polarized compact \kahler manifold of dimension $m$. We endow $L$ with a Hermitian metric $h$ with positive curvature. And we use $\omega=\frac{i}{2}\Theta_h$ as the \kahler form. By abuse of notation, we still use $h$ to denote the induced metric on the $k$-th power $L^k$. Then we have a Hermitian inner product on $H^0(X,L^k)$, defined by
$$<s_1,s_2>=\int_X h(s_1,s_2)\frac{\omega^m}{m!}$$
Let $\{s_i\}$ be an orthonormal basis of $\hcal_k=H^0(X,L^k)$. Then the  on-diagonal Bergman kernel, also called density of states function, is defined as
$$\rho_k(z)=\sum |s_i(z)|_h^2.$$

 The  on-diagonal Bergman kernel $\rho_k$ has very nice asymptotic expansion by
the results of Tian, Zelditch, Lu, etc. \cite{Tian1990On, Zelditch2000Szego, Lu2000On, Catlin, MM}. The asymptotic expansion has been found very useful, making $\rho_k$ an important bridge connecting \kahler geometry to algebraic geometry (see for example \cite{donaldson2001,DonaldsonSun2}).

Let $V\subset X$ be a closed submanifold of dimension $n$.
Let $\delta>0$ be a small number. We denote by $\hcal_k^{\delta k}\subset \hcal_k$ the subspace consisting of holomorphic sections whose vanishing orders along $V$ are no less than $\delta k$. The partial density function $\rho_k^\delta$ is then defined as the density of states function for $\hcal_k^{\delta k}$:
$$\rho_k^\delta(z)=\sum |s_i(z)|^2_h ,$$
	where $\{s_i\}$ is an orthonormal basis of $\hcal_k^{\delta k}$. When $V$ is a divisor,  the partial density functions have been studied by Zelditch-Zhou \cite{ZelditchZhou-Interface,ZelditchZhou2019}, Ross-Singer\cite{Ross2017} and Coman-Marinescu \cite{CM17} etc. It was first shown by Shiffman-Zelditch \cite{SZ2004} in the toric case, that there is a `forbidden region' in which the partial density function is exponentially small. Berman \cite{Berman2007} showed that there is an open forbidden region $R$ that asymptotically the partial density function is exponentially small on compact subsets and is equal to the usual density function on compact subsets of the complement of $\bar{R}$.
		In particular,  Ross-Singer\cite{Ross2017} studied the case when $V$ is smooth divisor, plus the condition that the line bundle and the metric is $S^1$-symmetric in a neighborhood $U$ of $V$. Let $\mu:U\to \R$ be the Hamiltonian of the $S^1$-action, normalized so that $\mu^{-1}(0)=V$. Then Ross-Singer proved a distributional asymptotic expansion for the partial density function and showed that the forbidden region is $\mu^{-1}[0,\delta)$:
	\begin{theo}[\cite{Ross2017}]
		For sufficiently small $\delta$, we have
		$$\rho_{k}^\delta \sim \left\{
		\begin{array}{lcl}
		O(k^{-\infty})& &\text{on } \mu^{-1}[0,\delta) \\
		\rho_k+O(k^{-\infty})& & \text{on } X\backslash \mu^{-1}[0,\delta]
		\end{array}
		\right.,$$
		where $\sim$ means  the equality holds on any given compact subset of $\mu^{-1}[0,\delta)$ and $X\backslash \mu^{-1}[0,\delta]$ respectively.
	\end{theo}
The term $k^{-\infty}$ in the theorem means a quantity depending on $k$ that is $O(k^{-m})$ as $k\rightarrow\infty$, for all $m\geq0$.
	We will adopt the notation used by Donaldson and denote $k^{-\infty}$ by $\epsilon (k)$. So, in order to avoid confusion, we use $\delta$ for the partial density function instead of $\epsilon$ used in \cite{Ross2017}.

\

In this article, we generalize Ross-Singer's result to the case when $V$ is a smooth subvariety. Let	$T^{m-n}=(S^1)^{m-n}$ denote the $(m-n)$-dimensional torus.
We assume that there is a neighborhood $U$ of $V$ that admits a holomorphic $T^{m-n}$-action. We also assume the line bundle $L$ and the metric $h$, hence $\omega$, are invariant under the $T^{m-n}$-action.

Let $\mu: U\to \R^{m-n}$ be the moment map of the torus action, normalized so that $\mu(V)=0$. We denote by $\mu_i$ the $i$-th component of $\mu$.
Then $\nu=\sum_{i=1}^{m-n}\mu_i:U\to \R$ is the moment map for the diagonal $S^1$-action.
Our main result is that the forbidden region is $\nu^{-1}([0,\delta))$:

\begin{theo}\label{main}
	For sufficiently small $\delta$, we have
	$$\rho_{k}^\delta \sim \left\{
	\begin{array}{lcl}
	\epk& &\text{on } \nu^{-1}[0,\delta) \\
	\rho_k+\epk& & \text{on } X\backslash \nu^{-1}[0,\delta]
	\end{array}
	\right.$$
\end{theo}
The proof uses the techniques we developed in \cite{SunSun,Annalen2019,Punctured}.
The idea of the proof is similar to that in \cite{Ross2017}. Unfortunately, we have not been able to give an asymptotic expansion for the renormalized partial Bergman kernel as Ross-Singer did in \cite{Ross2017}. It is also very interesting to know if there is still a forbidden region and what the boundary is, if the assumption of the torus-invariance is removed.

We would like to compare the partial density function to the logarithmic Bergman kernel studied by the author in \cite{sun2019logarithmic}, for which the vanishing order is $\geq 1$ instead of $\geq \delta k$. In \cite{sun2019logarithmic}, we showed that the asymptotics near the subvariety depend on the distance function to the subvariety. When $\delta$ is small enough, the moment function $\nu$ is comparable with the distance function squared. One can see the similarity between them.

 The asymptotics of off-diagonal Bergman kernel have been extensively used in the  value distribution theory of sections of line bundles by Shiffman-Zelditch and others\cite{bsz0,bsz1,bsz2,bsz3,sz4,sz1,sz2,sz3,SZ2004,szz,sz6,dsz1,dsz2,dsz3}. It would be very interesting and important to study the asymptotics of the partial off-diagonal Bergman kernel vanishing along a general subvariety.

	\

The structure of this article is as follows. We will first prove the key estimates on the complex plane. Then we apply the key lemma to the neighborhood $U$. Then with the help of the H\"{o}rmander's technique, we show that what we get on $U$ goes directly to $X$, hence proving the main theorem.

\

\textbf{Acknowledgements.} The author would also like to thank Professor Song Sun and Professor Chengjie Yu for many very helpful discussions.

	\section{One dimensional model}

Let $\omega=f\sqrt{-1}dz\wedge d\bar{z}$ be a smooth \kahler form on $\C$ invariant under the natural $S^1$-action. Let $\phi(|z|)$ be a smooth potential of $\omega$. By a linear change of coordinates, we can assume that $\phi=|z|^2+O(|z|^4)$. Let $\mu:\C\to\R$ be the moment map of the $S^1$-action normalized so that $\mu(0)=0$. We will use the notations $r=|z|$, $x=|z|^2$ and $t=-\log x$. So $\phi_{tt}=x\Delta \phi=xf$.
Then $$\mu=\int_0^x fdx.$$
And since $\lim_{t\to\infty}\phi_{t}=0$, we have
\begin{eqnarray*}
	\phi_t&=&-\int_{t}^{\infty}\phi_{tt}dt\\
	&=&\int_{x}^{0}fdx\\
		&=&-\mu
\end{eqnarray*}
Let $a\geq 0$ be an integer, $f_1(|z|)>0$ be a smooth function, we consider the $k$-th norm of the holomorphic function $z^a$:
$$I_a=\int_{\C}e^{-k\phi}|z|^{2a}f_1\sqrt{-1}dz\wedge d\bar{z}=2\pi\int_{-\infty}^{\infty}e^{-k\phi-(a+1)t+\log f_1}dt.$$
Denote by $g_a(t)=-k\phi-(a+1)t+\log f_1$. Then $g_a'=-k\phi_t-(a+1)+(\log f_1)_t$ and $g_a''=-k\phi_{tt}+(\log f_1)_{tt}$. Since $f_1=A+O(|z|^2)$, $A>0$, we have $(\log f_1)_{t}=O(\phi_{t})$ and $(\log f_1)_{tt}=O(\phi_{tt})$. So for $k$ large enough, $g_a$ is a concave function of $t$ which attains its only maximum at $t_a$ satisfying
$$\phi_t(t_a)=-\frac{a+1}{k}+\frac{(\log f_1)_t(t_a)}{k}.$$
Namely, $$\mu(t_a)=\frac{a+1}{k}+O(\frac{x_a}{k}),$$
when $x_a<1$.

\

We recall the following basic lemma used in \cite{SunSun}:

\begin{lem}\label{lemconcave}
	Let $f(x)$ be a concave function. Suppose $f'(x_0)<0$, then we have
	$$\int_{x_0}^\infty e^{f(x)}dx\leq\frac{e^{f(x_0)}}{-f'(x_0)}$$
\end{lem}

Let $0<R_1<R<1$ be two fixed numbers.
 When $\mu(t_a)>R$, we have $a>Rk-C$ for some $C$ independent of $k$, then $g_a'<0$ when $\mu\leq R$. Let $\mu(t_R)=R$ and $\mu(t_{R_1})=R_1$, then we let $t_m=\frac{t_R+t_{R_1}}{2}$. Clearly, $|g_a'(t_m)|>C_1k$ for some $C_1>0$ depending on $R$ and $ R_1$ but independent of $k$ and $a$. By the concavity of $g_a$, we have
 $$e^{g_a(t_m)-g_a(t_{R_1})}=\epk$$
 	Therefore, by lemma \ref{lemconcave}, we have
 	$$\int_{\mu<R_1}e^{g_a(t)}dt=\epk \int_{R_1\leq \mu\leq R}e^{g_a(t)}dt.$$

 	When $\mu(t_a)\leq R$, depending on $a$, there are two cases .

 	When $a\geq \sqrt{k}$, $g_a''(t_a)$ is large, so we can use Laplace's method to estimate the integral, namely
 	$$I_a\approx \frac{(2\pi)^{3/2}}{\sqrt{g_a''(t_a)}}.$$
 	More importantly, the  mass of the measure $e^{g_a(t)}dt$ is concentrated in a small neiborhood of $t_a$. More precisely, if $k\phi_{tt}(t_a)\delta_a^2>2\log^2 k$, then by lemma \ref{lemconcave}, we have $$I_a=(1+\epk)\int_{t_a-\delta_a}^{t_a+\delta_a}e^{g_a(t)}dt.$$
 	It is easy to see that there is a constant $C_2>0$ independent of $k$ such that if $\delta_a>C_2\frac{\log k}{\sqrt{a}}$, then the condition $k\phi_{tt}(t_a)\delta_a^2>2\log^2 k$ is satisfied.

	When $a<\sqrt{k}$, the term $k\phi_{tt}(t_a)$ may be too small to use Laplace's method. We have the following

	\begin{eqnarray*}
		\int_{x<k^{-1/4}} e^{g_a(t)}dt&\geq& e^{-(a+1)}e^{-(a+1)e}e^{g_a(t_a)}\\
		&\geq& e^{-(a+1)}e^{-(a+1)e}e^{-(a+1)+a\log\frac{a+1}{k}}\\
		&\geq& e^{-5(a+1)}(\frac{a+1}{k})^a
	\end{eqnarray*}
When $x=k^{-1/4}$, $\frac{e^{g_a(\log k^{1/4})}}{e^{-5(a+1)}(\frac{a+1}{k})^a}=\epk$. So by lemma
\ref{lemconcave}, we have for $a<\sqrt{k}$,
$$\int_{\mu\geq 2k^{-1/4}} e^{g_a(t)}dt	\leq \int_{x\geq k^{-1/4}} e^{g_a(t)}dt=\epk I_a$$

In summary, we have proved the following
\begin{lem}\label{1-dim-model}
For fixed $0<R_1<R<1$, there exist constants $C, C_1>0$ independent of $k$, such that for $k$ large enough we have the following estimations:
	\begin{itemize}
		\item[1)] If $a>Rk-C$, $$\int_{\mu<R_1}e^{g_a(t)}dt=\epk \int_{R_1\leq \mu\leq R}e^{g_a(t)}dt,$$
		where the $\epk$ term is independent of $a$.
			\item[2)] If $\sqrt{k}\leq a\leq Rk-C$,$$\int_{|\mu-\frac{a}{k}|\geq C_1\frac{\sqrt{a}\log k}{k}}e^{g_a(t)}dt=\epk \int_{|\mu-\frac{a}{k}|< C_1\frac{\sqrt{a}\log k}{k}}e^{g_a(t)}dt,$$
			where the $\epk$ term is independent of $a$.
				\item[3)] If $a<\sqrt{k}$, $$\int_{\mu\geq 2k^{-1/4}} e^{g_a(t)}dt	=\epk \int_{\mu< 2k^{-1/4}} e^{g_a(t)}dt,$$
				where the $\epk$ term is independent of $a$.
	\end{itemize}

\end{lem}

	\section{proof of the main theorem}
By assumption, there is a neighborhood $U$ of $V$ that is covered with charts that admit standard coordinates $(z,w)$ so that $V=\{z=0\}$, where $z=(z_1,\cdots,z_{m-n})\in \C^{m-n}$ and $w=(w_1,\cdots,w_n)\in \C^n$, and
$$\lambda\cdot(z,w)=(\lambda z, w),$$
for $\lambda=(\lambda_1,\cdots,\lambda_{m-n})$, where $\lambda z=(\lambda_1z_1,\cdots,\lambda_{m-n}z_{m-n})$. We have the following lemma,
\begin{lem}
	There is a $T^{m-n}$-invariant biholomorphism $\Psi$ from the neighborhood $U$ to an open neighborhood of $V$ in the total space of the vector bundle $L_1\oplus L_2\cdots \oplus L_{m-n}$, where each $L_i$ is a holomorphic line bundle.
\end{lem}
 The proof is very similar to that of lemma 5.3 in \cite{RossSinger}, and for the convenience of the reader, we include a proof here.
\begin{proof}
	Let $(z_{\alpha},w_{\alpha})$ and $(z_\beta,w_\beta)$ be two sets of standard coordinates. Then the transition functions between them are necessarily of the form
	$$z_\beta=\lambda_{\alpha\beta}(z_\alpha,w_\alpha)z_\alpha,\quad w_\beta=\tau_{\alpha\beta}(z_\alpha, w_\alpha) $$
	where $\lambda_{\alpha\beta}$ and $\tau_{\alpha\beta}$ are holomorphic and take values in $\C^{*(m-n)}$ and $\C^{(m-n)}$ respectively.
	We understand the torus $T^{m-n}$-action on the $w$-part is trivial. Then both $\lambda_{\alpha\beta}$ and $\tau_{\alpha\beta}$ are $T^{m-n}$-invariant, so their dependence on $z_\alpha$ are trivial. So the transition functions
	$$z_\beta=\lambda_{\alpha\beta}(w_\alpha)z_\alpha$$
	defines a holomorphic vector bundle $E$ of rank $m-n$ over $V$, with structure group $T^{m-n}$. Therefore, $E$ decomposes as $E=L_1\oplus L_2\cdots\oplus L_{m-n}$, and $U$ is biholomorphic to an open neighborhood of $V$ in $E$. Clearly, the biholomorphism is $T^{m-n}$-invariant.
\end{proof}
We consider $U$ as an open set in $\oplus L_i$ and denote by $\pi:U\to V$ the projection. Then for $A$ small enough, $U$ contains the poly-disc bundle $U_A=\cap_{1\leq i\leq m-n}\mu_i^{-1}(x<A)$. For each $p\in V$, the fiber of $U_A$ is an open poly-disc $D_A$ in $\C^{m-n}$. By the torus-invariance, the restriction $\phi_p$ of $\phi$ depends only on $|z_1|,\cdots,|z_{m-n}|$, so does the restriction $\omega_p$ of $\omega$. So the monomials $z^a=\prod z_i^{a_i}$ are all orthogonal to each other on $D_p$, where  $a=(a_1,\cdots,a_{m-n})$ is an $(m-n)$-tuple of non-negative integers.

\begin{lem}
	Let $W_\delta=\{q\in D_p|\nu(q)<\delta-k^{1/4} \}$.
	When $\delta$ is small enough, for each holomorphic function $f=\sum_{\sum a_i\geq \delta k}c_az^a$, we have
	$$\int_{W_\delta}|f|^2e^{-k\phi_p}\omega_p^{m-n}=\epk\int_{D_p}|f|^2e^{-k\phi_p}\omega_p^{m-n}$$
\end{lem}
\begin{proof}
	By orthogonality, $\int |f|^2e^{-k\phi_p}\omega_p^{m-n}=\sum_{\sum a_i\geq \delta k} |c_a|^2\int |z^a|^2 e^{-k\phi_p}\omega_p^{m-n}$. So it suffices to prove the lemma for the monomials $z^a$ satisfying $\sum a_i\geq \delta k$.

	By making $A$ a little smaller, a compactness argument shows that the constants $C$ and $C_1$ in lemma \ref{1-dim-model} can be chosen to be the same for all $z_i$-discs. For each point $p$ in the open set , at least one $1\leq i\leq m-n$ satisfies the condition that $\mu_i< \frac{a_i}{k}-k^{-3/8}$. Therefore, lemma \ref{1-dim-model} shows that the mass of $z^a$ in $W$ is $\epk$ relative to the mass in $D_A\backslash W_\delta$. And the theorem is proved

\end{proof}
As a corollary, we can prove the part about the forbidden region of theorem \ref{main}.
\begin{proof}[Proof of the first half of theorem \ref{main}]
	We can take an arbitrary orthonormal basis $\{s_i\}$ for the space $\hcal^{\delta k}_k$. Then the lemma above shows that each $s_i$ has $L^2$-mass of the size $\epk$ in the region $W_\delta$. We reall that for each point $p\in X$, there is an unit section $s_p$, called the peak section, satisfying $|s_p(p)|^2_h=\rho_k(p)$ and $s_p\perp s$ as long as $s(p)=0$.

We take the inner product of $s_i$ with the peak section$s_p$ for each $p\in W$
	$$|s_i(p)|=|<s_i,s_p>|\sqrt{\rho_k(p)}\leq \parallel s_i\parallel \sqrt{\rho_k(p)} $$
	 we see that the point-wise norm $|s_i(p)|$ is also $\epk$.   Since the dimension of $\hcal^{\delta k}_k$ is $O(k^m)$. We get the conclusion.
\end{proof}

We denote by $T_i$ the $i$-th component of $T^{m-n}$.
For each $L_i$, we fix a $T_i$-invariant section $s_i$ such that $\mu_i^{-1}(0)=\{s_i=0\}$. So we can choose local frame $e_i$ of $L_i$, so that $s_i=z_i e_i$ in the standard torus-invariance coordinates.

For each $s\in H^0(V,kL-\otimes L_i^{a_i})$, where $a_i\geq 0$, $\pi^*s\in H^0(U,kL-\otimes L_i^{a_i})$, so $\tilde{s}=\pi^*s\otimes\prod s_i^{a_i}\in H^0(U,L^k)$.
 We denote by $G_\delta$ the span of all such sections for $\sum_{i=1}^{m-n}a_i<\delta k$. Then clearly, for each section $\beta\in \hcal^{\delta k}_k$, we have $\beta|_U\perp G_\delta$. Moreover, since each $\tilde{s}$ is represented by a holomorphic function of the form
 $$\sum_{\sum_{i=1}^{m-n}a_i<\delta k} f_a(w)z^a,$$
  by lemma \ref{1-dim-model}, each section $\tilde{s}\in G_\delta$ has its mass concentrated within  the set $Y_\delta=\{q|\nu(q)\leq \delta+2(m-n)k^{-1/4} \}$ and decays fast as $\nu$ gets bigger. More precisely, the measure outside $Y_\delta$ is $\epk$ relative to that inside $Y_\delta$. To connect $G_\delta$ to $(\hcal^{\delta k}_k)^\perp$, we need to use H\"ormander's $L^2$ technique. The following lemma is well-known, see for example \cite{Tian1990On}.

\begin{lem}
	Suppose $(M,g)$ is a complete \kahler manifold of complex dimension $n$, $\mathcal L$ is a line bundle on $M$ with hermitian metric $h$. If
	$$\langle-2\pi i \Theta_h+Ric(g),v\wedge \bar{v}\rangle_g\geq C|v|^2_g$$
	for any tangent vector $v$ of type $(1,0)$ at any point of $M$, where $C>0$ is a constant and $\Theta_h$ is the curvature form of $h$. Then for any smooth $\mathcal L$-valued $(0,1)$-form $\alpha$ on $M$ with $\bar{\partial}\alpha=0$ and $\int_M|\alpha|^2dV_g$ finite, there exists a smooth $\mathcal L$-valued function $\beta$ on $M$ such that $\bar{\partial}\beta=\alpha$ and $$\int_M |\beta|^2dV_g\leq \frac{1}{C}|\alpha|^2dV_g$$
	where $dV_g$ is the volume form of $g$ and the norms are induced by $h$ and $g$.
\end{lem}

Now let $\chi\in C^{\infty}_0(U)$ be a non-negative cut-off function, that equals 1 in a neighborhood of the subset $\{p\in U|\nu(p)\leq \delta+c \}$ for some small $c>0$. Then for each $\tilde{s}\in G_\delta$ with $\parallel \tilde{s}\parallel^2=1$, we take the orthogonal projection of $\chi  \tilde{s}$ onto $\hcal_k$. The way to do so is to solve the $\dbar$-equation
$$\dbar v=\dbar \chi \otimes \tilde{s}.$$
And since $|\tilde{s}|$ is $\epk$ on the place where $\dbar\chi\neq 0$, we can apply the lemma above to see that the global holomorphic section $\chi\tilde{s}-v\in \hcal_k$ has $\epk$ $L^2$-mass outside $U$ and the $L^2$-mass of $\chi\tilde{s}-v-\tilde{s}$ within $U$ is also $\epk$.
So for each $f\in \hcal^{\delta k}_k$ with $\parallel f\parallel^2=1$, we have
$$|<f,\chi\tilde{s}-v>|=\epk $$
 We then take the orthogonal projection of $\chi\tilde{s}-v$ onto $(\hcal^{\delta k}_k)^\perp$ and denote by $\bar{s}$ the section obtained. So we have
 $$\int_U |\bar{s}-\tilde{s}|^2e^{-k\phi}\omega^m=\epk.$$
We denote by $J$ the map we just constructed from $G_\delta$ to $(\hcal^{\delta k}_k)^\perp$. So $J$ is almost an isometry, namely
$$\parallel J(\tilde{s})\parallel^2=(1+\epk)\parallel \tilde{s}\parallel^2.$$
Therefore, if we start with an orthonormal basis of $G_\delta$, we obtain an almost orthonormal basis of $(\hcal^{\delta k}_k)^\perp$. And then after a Gram-Schmidt process, we get an orthonormal basis. We denote by $\rho_G$ the density function of $G_\delta$, considered as a function on $X$ by extenstion by zero. Since the density function of $(\hcal^{\delta k}_k)^\perp$ is $\rho_k-\rho_k^{\delta k}$, we have
$$\int_X (\rho_k-\rho_k^{\delta }-\rho_G)\omega^m=\epk.$$
Since, $\rho_G$ has $\epk$ mass outside  $Y_\delta$, we have
$$\int_{X\backslash Y_\delta} (\rho_k-\rho_k^{\delta })\omega^m=\epk.$$
Recall that by a local construction of the peak section $s_p$ for any $p\in X$(see for example\cite{Donaldson15}), $|s_p|$ decays exponentially away from $p$. More precisely, $|s_p(q)|=\epk$, when $d(p,q)>\frac{\log k}{\sqrt{k}}$. When $\delta$ is small enough, $\nu$ is comparable with the distance function $d^2(p,q)$, so if we denote by
 $Y'_\delta=\{q|\nu(q)\leq \delta+3(m-n)k^{-1/4} \}$, then for any unit $s\in (\hcal^{\delta k}_k)^\perp$, we have
 $$|<s,s_p>|=\epk,$$
for $p\in X\backslash Y'_\delta$. Therefore $|s(p)|=\epk$, $p\in X\backslash Y'_\delta$. Since $\lim_{k\to \infty }Y'_\delta=\{q|\nu(q)\leq \delta \}$, this implies the second part of theorem \ref{main}.

	\bibliographystyle{plain}


\section{Introduction}
Recent measurements by ALICE~\cite{ALICE:2017jyt} show a smooth evolution of strange to non-strange particle ratios across colliding systems (pp, p--Pb, and Pb--Pb) as a function of charged-particle multiplicities. The observed enhancement scales with the strangeness content of the hadron. Similar features of various particle ratios measured across different collision systems are also observed in Ref.~\cite{Acharya:2018orn}. Collective-like effects have also been observed in small systems (pp and p-Pb collisions), but there is so far no evidence for jet quenching in small systems~\cite{Nagle:2018nvi}. In popular event generators such as PYTHIA8~\cite{Sjostrand:2014zea} the observed effects are only reproduced when novel final-state interactions (ropes and shoving) are introduced. In EPOS-LHC~\cite{Pierog:2013ria}, one can qualitatively describe the effect as by inclusion of QGP production in small systems. Small systems in EPOS-LHC are a mix of corona (QCD) and core (QGP), and it is the relative fraction of these that mainly drives the change with system size. To understand the origins of these phenomena, event shape observables such as the unweighted transverse spherocity ($S_{\rm 0}^{p_{\rm T = 1}}$) and the relative transverse activity classifier ($R_{\rm{T}}$) can be exploited as a powerful tools to disentangle soft and hard particle production.

In this contribution, we report the production of light-flavor hadrons for various $S_{\rm 0}^{p_{\rm T} = 1}$ classes in pp collisions at $\sqrt{s}$ = 13 $\textrm{TeV}$. In addition, the system size dependence of charged particle production in pp, p--Pb, and Pb--Pb collisions at  $\sqrt{s_{\rm NN}}$ = 5.02 TeV as a function of $R_{\rm{T}}$ is reported. The evolution of $\langle p_{\rm T}\rangle$ in different topological regions as a function of $R_{\rm{T}}$ are presented. Finally, we present a search for jet quenching behavior in small collision systems using similar techniques.

\section{Event shape observables}
Event shape observables have the capability to separate events with back-to-back jet structures and events dominated by multiple soft scatterings. Here, we use $S_{\rm 0}^{p_{\rm T} = 1}$ and $R_{\rm T}$ as estimators to disentangle events with these topological limits.

\subsection{Unweighted transverse spherocity ($S_{\rm 0}^{p_{\rm T} = 1}$)}

The unweighted transverse spherocity, inspired from transverse spherocity~\cite{Acharya:2019mzb}, is given by,
\begin{eqnarray}
S_{0}^{p_{\rm T} = 1} = \frac{\pi^{2}}{4} \min_{\hat{n}} \bigg(\frac{\Sigma_{i}~|p_{\rm T_{i}}\times\hat{n}|}{N_{\rm trks}}\bigg)^{2},
\label{eq1}
\end{eqnarray}
where $N_{\rm trks}$ is total number of charged-particle tracks in a given event and $\hat{n}$ is a unit vector that minimizes Eq.~\ref{eq1}. Here, $S_{0}^{p_{\rm T} = 1}$ is calculated using charged-particle tracks that have $p_{\rm T} > $ 0.15 GeV/$c$ for events with at least 10 charged-particle tracks to ensure that the concept of a topology is statistically meaningful. The charged-particle tracks are reconstructed using the Time Projection Chamber (TPC), within the pseudorapidity interval $|\eta|<0.8$. Unlike the estimator discussed in Ref.~\cite{Acharya:2019mzb}, the $p_{\rm T}$ of each track is normalized to unity ($p_{\rm T}$ = 1). This modification in the estimator is required to minimize biases which affect neutral particle yields. By definition, the $S_{0}^{p_{\rm T = 1}}$ estimator varies between the values 0 and 1, where the two extreme limits correspond to the two different topological limits. Events with $S_{0}^{p_{\rm T} = 1} \rightarrow$ 0 are dominated by a single back-to-back jet-like structure while events with $S_{0}^{p_{\rm T} = 1} \rightarrow$ 1 are dominated by isotropic particle production. From here onwards, the events within the bottom 20\% quantile of the $S_{0}^{p_{\rm T} = 1}$ distribution are referred to as the jetty events while the top 20\% quantile of the $S_{0}^{p_{\rm T} = 1}$ distribution are referred to as the isotropic events. As reported by ALICE~\cite{ALICE:2017jyt}, the strangeness enhancement is observed in high-multiplicity pp collisions, in addition to the $S_{0}^{p_{\rm T} = 1}$ selection, a top 10\% high-multiplicity requirement is also imposed for the event selection (multiplicity class I-III~\cite{Acharya:2019mzb}). The multiplicity is estimated using two sub-detector systems in ALICE. The multiplicity is estimated at mid-rapidity ($|\eta|<0.8$) using the Inner Tracking System (ITS) and it is referred as $N_{\rm SPD}$ estimator. The multiplicity is also estimated using the V0M scintillators which measure the charged-particle multiplicity at forward rapidity (2.8 $<\eta<$ 5.1 and -3.7 $<\eta<$ -1.7) and  it is referred as V0M estimator. The major difference in the measurements of particle production using two different estimators is that for V0M-based selections, the pseudo-rapidity regions for multiplicity estimation and for the measurement of the particle spectra are different, while for $N_{\rm SPD}$-based selections, the multiplicity and particle spectra are measured in the same rapidity region.

\subsection{Relative transverse activity classifier ($R_{\rm T}$)}
The measurements of particle production as a function of relative transverse activity classifier allow for the estimation of the production of hard probes relative to the underlying event. To ensure that a hard scattering took place in the event, analysed events are required to have a leading trigger particle of 8 $< p_{\rm T} <$ 15 GeV/$c$. An event is classified into three different azimuthal regions, relative to the leading track. Assuming $\phi_{\rm trig.}$ as the azimuthal angle for the leading trigger particle and $\phi_{\rm assoc.}$ as the azimuthal angle of the associated particles, the three regions are classified as the following,
\begin{itemize}
\item Near side: $|\phi_{\rm trig.} - \phi_{\rm assoc.}| < \frac{\pi}{3}$
\item Away side: $|\phi_{\rm trig.} - \phi_{\rm assoc.}| > \frac{2\pi}{3}$
\item Transverse side: $\frac{\pi}{3} \leq |\phi_{\rm trig.} - \phi_{\rm assoc.}| \leq \frac{2\pi}{3}$
\end{itemize}
The particle production in near side is expected to be dominated by jet fragmentation while the away side region captures some of the back-scattered jets. Both regions are also affected by the underlying event. The transverse region is mostly dominated by the underlying event, with possible contamination from initial- and final-state radiation. This leading-$p_{\rm T}$ selection of  8 $< p_{\rm T} <$ 15 GeV/$c$ ensures that the number density in the transverse region remains nearly independent of leading particle $p_{\rm T}$~\cite{Acharya:2019nqn}. This requirement also ensures that the effects of elliptic flow are negligible, which is important as the goal of the approach is to search for jet quenching effects in small collision systems. The relative transverse activity classifier ($R_{\rm T}$) is defined as,
\begin{eqnarray}
R_{\rm T} = \frac{N_{\rm ch}^{\rm TS}}{\langle N_{\rm ch}^{\rm TS} \rangle},
\label{eq2}
\end{eqnarray}
where $N_{\rm ch}^{\rm TS}$ is the charged-particle multiplicity in the transverse region~\cite{Martin:2016igp,Ortiz:2017jaz}. Using $R_{\rm T}$, one can vary the magnitude of underlying event. For example, as $R_{\rm T} \rightarrow$ 0 are the events are populated with less underlying events and they are expected to be compatible to the jetty events as defined from $S_{\rm 0}^{p_{\rm T} = 1}$.

\section{Results and Discussion}

\subsection{Identified particle production as a function of $S_{\rm 0}^{p_{\rm T}= 1}$}

\begin{figure}[h]
\centering
\includegraphics[width=14pc]{2019-10-31-2019-10-31-V0M_Pion_Spectra.pdf}
\includegraphics[width=14pc]{2019-10-31-2019-10-31-Trks_Pion_Spectra.pdf}
\caption{\label{fig1} Transverse momentum spectra of pions for different $S_{\rm 0}^{p_{\rm T}= 1}$ classes in pp collisions at $\sqrt{s}$ = 13 TeV. The left (right) panel shows the $p_{\rm T}$ spectra where the multiplicity selection was done using the V0M ($N_{\rm SPD}$) multiplicity estimator. The bottom panels show the ratio of $p_{\rm T}$-spectra from isotropic and jetty events to the $S_{\rm 0}^{p_{\rm T}= 1}$-integrated events.}
\end{figure}

Figure~\ref{fig1} shows the transverse momentum spectra of pions for different $S_{\rm 0}^{p_{\rm T}= 1}$ classes in pp collisions at $\sqrt{s}$ = 13 TeV. The left (right) panel shows the $p_{\rm T}$ spectrum where the multiplicity selection was done using the V0M ($N_{\rm SPD}$) multiplicity estimator. Pions are identified from the specific energy loss ($\rm{d}E/\rm{d}x$) measured by the TPC and time-of-flight (TOF) using the TOF detector. The bottom panels show the ratio of $p_{\rm T}$-spectra from events identified as isotropic and jetty to the $S_{\rm 0}^{p_{\rm T}= 1}$-integrated events. A clear dependence of the $p_{\rm T}$-spectra as a function of $S_{\rm 0}^{p_{\rm T} = 1}$ is observed in the results from both the estimators. However, as can be seen in Fig.~\ref{fig2} (pion $\langle p_{\rm T}\rangle$ vs integrated yield), the V0M-triggered isotropic and jetty events have larger multiplicity difference while the $N_{\rm SPD}$-triggered events disentangle soft and hard events more accurately in small multiplicity gap.
\begin{figure}[h]
\centering
\includegraphics[width=15pc]{2019-10-31-dndy.pdf}
\begin{minipage}[b]{14pc}\caption{\label{fig2} Pion $\langle p_{\rm T}\rangle$ vs integrated yield for different multiplicity estimators in different $S_{\rm 0}^{p_{\rm T}= 1}$ classes for pp collisions at $\sqrt{s}$ = 13 TeV.}
\end{minipage}
\end{figure}
Figure~\ref{fig3} shows proton-to-pion (left) and $\Xi$-to-pion (right) ratio as a function of $p_{\rm T}$ in different $S_{\rm 0}^{p_{\rm T} = 1}$ classes for pp collisions at $\sqrt{s}$ = 13 TeV, where the multiplicity selection was done using the $N_{\rm SPD}$ multiplicity estimator. The bottom panels show the double ratio of particle ratios from isotropic and jetty events to the $S_{\rm 0}^{p_{\rm T} = 1}$-integrated events. Here, protons are identified using TPC and TOF detectors. The $\Xi$ yields are extracted from invariant mass distributions of its identified decay products after topological selections. The proton-to-pion ratio show a larger enhancement in the intermediate-$p_{\rm T}$ region for isotropic events compared to jetty events, which seems to be reminiscent of a similar effect observed in Pb-Pb collisions~\cite{Acharya:2018orn}. The $\Xi$-to-pion ratio suggests that the strange particle production is higher in isotropic events compared to jetty events. In general, the particle ratios are not described by PYTHIA8 and EPOS-LHC. However, the trend of the double ratios are better described by both PYTHIA8 and EPOS-LHC.
\begin{figure}[h]
\centering
\includegraphics[width=18pc]{2019-10-31-Trks_Proton_to_Pion_ratio_MC_Predictions.pdf}
\includegraphics[width=18pc]{2019-10-31-Trks_Xi_to_Pion_ratio_MC_Predictions.pdf}
\caption{\label{fig3} Proton-to-pion (left) and $\Xi$-to-pion (right) ratio as a function $p_{\rm T}$ in different $S_{\rm 0}^{p_{\rm T}= 1}$ classes for pp collisions at $\sqrt{s}$ = 13 TeV where the multiplicity selection was done using the $N_{\rm SPD}$ multiplicity estimator. The bottom panels show the double ratio of particle ratios from isotropic and jetty events to the $S_{\rm 0}^{p_{\rm T}= 1}$-integrated events. The ratios are compared with predictions from PYTHIA 8 and EPOS-LHC.}
\end{figure}

\subsection{System size dependence of charge particle production as a function of $R_{\rm T}$}
Figure~\ref{fig4} shows the system size dependence of $\langle p_{\rm T}\rangle$ of charged-particles as a function of $R_{\rm T}$ in the near (left), away (middle), and transverse (right) regions. The near and away side $\langle p_{\rm T}\rangle$ for pp and p--Pb collisions decreases at low-$R_{\rm T}$ and saturates for high-$R_{\rm T}$. The $\langle p_{\rm T}\rangle$ for transverse side increases with $R_{\rm T}$. The contribution from the near and away side jet dominates at low-$R_{\rm T}$ and the values are similar for all systems as one would naively expect for $R_{\rm T}\rightarrow$ 0. For large $R_{\rm T}$, the $\langle p_{\rm T}\rangle$ approaches a similar value in all three topological regions for a given system.

\begin{figure}[h]
\centering
\includegraphics[width=35pc]{2020-05-22-Figure2a_0.pdf}
\caption{\label{fig4} System size dependence of $\langle p_{\rm T}\rangle$ of charged-particles as a function of $R_{\rm T}$ in the near (left), away (middle), and transverse (right) regions.}
\end{figure}
To explore the presence of jet-quenching effects we have calculated the $I_{\rm pp,p-Pb,Pb-Pb}$, an observable which is calculated from the yields of different topological regions, as a function of $\langle N_{\rm ch}^{\rm TS} \rangle$ for different V0M/V0A multiplicity classes of pp, p--Pb, and Pb--Pb collisions. The $I_{\rm pp,p-Pb,Pb-Pb}$ is an analogous quantity calculated as in Ref.~\cite{Aamodt:2011vg}. The $I_{\rm pp,p-Pb,Pb-Pb}$ is sensitive to medium effects and any suppression in the away side would indicate the presence of jet quenching, while an enhancement in the near side would indicate the presence of medium effects and a bias due to trigger particle selection.  It is defined as the ratio of yield in the near or away region (after subtraction of underlying events) in different collision systems to the yield in the near or away region in minimum bias pp collisions. It can be expressed as,

\begin{equation}
I_{\rm pp,p-Pb,Pb-Pb} = \frac{Y^{\rm pp,p-Pb,Pb-Pb} - Y^{\rm pp,p-Pb,Pb-Pb}_{\rm TS}}{Y^{\rm pp~min. bias} - Y^{\rm pp~min. bias}_{\rm TS}}.
\end{equation}
Here, $Y$ represents the integrated yield of charged particles. Note that, for these results there is no direct selection on $N_{\rm ch}^{\rm TS}$, as direct selection on $N_{\rm ch}^{\rm TS}$ biases the near and away side yields~\cite{Ortiz:2020dph}. Thus, the events are selected based on the forward rapidity estimator (V0M for pp and Pb--Pb collisions and V0A for p--Pb collisions) and the corresponding $N_{\rm ch}^{\rm TS}$ are calculated. Figure~\ref{fig5} shows the $I_{\rm pp,p-Pb,Pb-Pb}$ for the range 4 $< p_{\rm T}^{assoc.} < $ 6 GeV/$c$ as a function of $\langle N_{\rm ch}^{\rm TS} \rangle$ in different V0M/V0A multiplicity classes for the near (left) and away (right) side in pp, p--Pb, and Pb--Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV. The values of $I_{\rm PbPb}$ for central and peripheral Pb-Pb collisions show a similar trend for both the near and away sides as reported by ALICE in Ref.~\cite{Aamodt:2011vg} at Pb--Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV. For small collision systems, no enhancement (suppression) of $I_{\rm pp,p-Pb}$ is observed in near (away) sides for pp or p--Pb collisions. This may indicate the absence of jet-quenching effects in small collision systems for the measured  $\langle N_{\rm ch}^{\rm TS} \rangle$ ranges.

\begin{figure}[h]
\centering
\includegraphics[width=16pc]{2020-05-22-Figure4_IAA_NS.png}
\includegraphics[width=16pc]{2020-05-22-Figure4_IAA_AS.png}
\caption{\label{fig5} $I_{\rm pp,p-Pb,Pb-Pb}$ as a function of $\langle N_{\rm ch}^{\rm TS} \rangle$ in different V0M/V0A multiplicity classes for the near (left) and away (right) regions in pp, p--Pb, and Pb--Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV.}
\end{figure}

\section{Summary}
In summary, using the event shape observables $S_{\rm 0}^{p_{\rm T} = 1}$ and $R_{\rm T}$, the magnitude of the underlying events can be varied to study the events separately with different topological limits (jetty vs isotropic). A clear dependence of light flavor particle $p_{\rm T}$ spectra on $S_{\rm 0}^{p_{\rm T}= 1}$ is observed. The system size study of charged-particle production suggests that the contribution from the near and away side jet dominates at low-$R_{\rm T}$ and the values are similar for all systems, as one would naively expect for $R_{\rm T}\rightarrow$0. For large $R_{\rm T}$, the $\langle p_{\rm T}\rangle$ approaches a similar value in all three topological regions for a given system. In contrast to Pb--Pb collisions, no enhancement (suppression) of $I_{\rm pp,p-Pb}$ is observed in the near (away) sides for pp, p--Pb collisions, which indicates the absence of jet-quenching effects in small collision systems for the measured  $\langle N_{\rm ch}^{\rm TS} \rangle$ ranges.

\section*{Acknowledgments}
S.T. acknowledges the support from CONACyT under the Grant No. A1-S-22917 and postdoctoral fellowship of DGAPA UNAM.

\section*{References}


\section{Introduction}

Differential equations are prevalent in every field of science and engineering, ranging from physics to economics. Thus, extensive research has been done on developing numerical methods for solving differential equations.
With the unprecedented availability of computational power, neural networks hold promise in redefining how computational problems are solved. Among other applications, unsupervised neural networks are capable of  efficiently   solving  differential equations \cite{ourLib2020, pnas2018,  lagaris1998, nips2018,  mattheakis2020hamiltonian}.
These are unsupervised, data-free methods where the  optimization depends solely on the network predictions.
The neural network solvers pose several advantages over numerical integrators: the obtained solutions are analytical and differentiable \cite{lagaris1998},
{networks are more robust to} the `curse of dimensionality' \cite{pnas2018},  numerical errors are not accumulated \cite{mattheakis2020hamiltonian}, and a family of solutions corresponding to different initial or boundary conditions can be constructed \cite{cedric2020}.

Differential eigenvalue equations with boundary conditions
appear in a wide range of problems, including quantum mechanics and applied mathematics.
Efficient numerical iterative methods, such as finite difference method,  have been developed  for solving eigenvalue problems, but they share the drawbacks common to all numerical integrators.
{Lagaris et al.} \cite{lagaris_eigen} have shown that neural networks are able to solve eigenvalue problems and proposed a  {  partially iterative method that solves a differential equation with a fixed eigenvalue at each iteration.}
{Our contribution includes} a novel unsupervised neural network architecture that  simultaneously learns eigenvalues and the associated eigenfunctions using a scanning mechanism.
The proposed technique is an extension to neural network differential equation solvers and, consequently, acquires all the benefits that  network solvers have over  numerical integrators. Moreover, our method has an additional advantage over integrators in that it discovers solutions that identically satisfy the boundary conditions.
We assess the performance of the proposed architecture by solving two  standard eigenvalue problems of quantum mechanics, namely, the infinite well, and the quantum harmonic oscillator.

\section{Methodology}
\label{gen_inst}

We consider an eigenvalue problem that exhibits the form:
\begin{align}
\label{eq:EigenProblem}
    \mathcal{L}f(x) = \lambda f(x),
\end{align}
where $x$ is the spatial variable, $\mathcal{L}$ is a differential operator which depends on $x$ and its derivatives, $f(x)$ is the eigenfunction, and $\lambda$ is the associated eigenvalue. We assume homogeneous  Dirichlet boundary conditions at $x_{L}$ and $x_{R}$ such that $f(\xL)=f(\xR)=f_b$,
where $f_b$ is a constant boundary value.
For a given and fixed eigenvalue, Eq. (\ref{eq:EigenProblem}) is an  equation that  can be solved by neural network methods suggested in \cite{ourLib2020, lagaris_eigen,lagaris1998}.
We introduce a new  architecture shown in Fig. \ref{fig:finalarch}, which is capable of solving Eq. (\ref{eq:EigenProblem}) when both $f(x)$ and $\lambda$ are unknown.
The network  takes two inputs,  $x$ and  $1$. The constant input feeds a single linear neuron (affine transformation) that is
updated through optimization, allowing  the network to find constant   $\hlam$.
The $x$ and $\hlam$ feed a  feed-forward fully-connected network that returns an output function $N(x,\hlam)$. The predicted eigenfunctions $f(x,\hlam)$ is defined by using  a parametric trick, similar to Ref. \cite{mattheakis2020hamiltonian}, according to the equation:
\begin{align}
\label{eq:parSol}
    f(x,\hlam) = f_b + g(x)N(x,\hlam),
\end{align}
where we employ the parametric function
\begin{align}
\label{eq:parametric}
    g(x) = \left(1-e^{-(x-\xL)}\right)\left(1-e^{-(x-\xR)}\right),
\end{align}
which embeds the boundary conditions in the $f(x,\hlam)$.
\begin{figure}[h]
    \centering
    \includegraphics[width = 0.7\textwidth]{/arch3final.pdf}
    \caption{Adopted eigenvalue problem architecture. }
    \label{fig:finalarch}
\end{figure}

Our  aim   is to discover pairs of $f(x, \hlam )$ and $\hlam$  that satisfy Eq. (\ref{eq:EigenProblem}). This is achieved by minimizing, during the optimization, a loss function $L$   defined by Eq. (\ref{eq:EigenProblem}) as:
\begin{align}
\label{eq:Loss}
    L &=  L_\text{DE} +  L_\text{reg} \nonumber \\
      &=  \Big \langle \big( \mathcal{L}f(x, \lambda) - \lambda f(x, \lambda) \big)^2  \Big \rangle_x +  L_\text{reg},
\end{align}
 where $\langle \cdot \rangle_x$ represents  averaging with respect to $x$. Any derivative with respect to  $x$ contained in  $\mathcal{L}$ is calculated by using the auto-differentiation technique \cite{Paszke2017AutomaticDI}.
The $L_\text{reg}$  in Eq. (\ref{eq:Loss}) contains regularization loss functions and is defined as: $L_\text{reg} = \nu_f L_f +\nu_\lambda L_\lambda + \nu_\text{drive} L_\text{drive}$. { Empirically, for the problems discussed below, we found  the optimal set  $\nu_ f= \nu_\lambda=\nu_\text{drive}=1$.} The $L_f$ and $L_\lambda$  are used to avoid learning  trivial eigenfunctions and eigenvalues respectively, while $ L_\text{drive}$  motivates the network to scan for higher eigenvalues, {as we explain below}. The regularization functions   are defined as:
\begin{align}
\label{eq:nontriv}
    L_{f} = \frac{1}{f(x, \lambda)^2}, \quad \quad \quad \quad
    L_{\lambda} = \frac{1}{\lambda^2}, \quad \quad \quad \quad
    L_{\text{drive}} = e^{-\lambda + c}.
\end{align}
During the optimization, a scheduled scanning algorithm increases  the $L_\text{drive}$ by increasing $c$  in regular intervals. That forces the network to search for larger eigenvalues and the associated eigenfunctions.  At each interval the network is optimized and  the model parameters are stored when sufficiently low $L_\text{DE}$ is achieved.
We emphasize that the loss function solely  depends  on the predictions of the network and, therefore, the training process is data-free, resulting in an unsupervised learning method.
For the training, a batch of $x$ points in the interval $[\xL, \xR]$ is selected as input. In every training iteration (epoch) the input points are perturbed by a Gaussian noise \cite{mattheakis2020hamiltonian}. Adam optimizer is used \cite{adam} with a learning rate of $8 \cdot 10^{-3}$.
We use two hidden layers of 50 neurons per layer with  trigonometric $\sin(\cdot)$  activation function.  The use of  $\sin(\cdot)$ instead of more common activation functions, such as Sigmoid$(\cdot)$ and $\tanh(\cdot)$, significantly accelerates the network's convergence to a solution \cite{mattheakis2020hamiltonian}.  We  implemented  the proposed neural network in pytorch  \cite{Paszke2017AutomaticDI} and published the code on github \footnote{https://github.com/henry1jin/eigeNN}.

\section{Experiments}
\label{headings}

We evaluate the effectiveness of the proposed method  by solving the eigenvalue problem defined by Schrodinger's equation. This is a fundamental equation in quantum mechanics that  describes the state wavefunction $\psi(x)$  and energy $E$ of a quantum system. We are interested in solving the one-dimensional stationary Schrodinger equation  defined as:
\begin{align}
\label{eq:QMeigen}

    \left[ -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x)\right] \psi(x) &= E\psi(x),
\end{align}
where $\hbar$ and $m$ stand for Planck constant and mass {which}, without loss of  generality, can be set to   $\hbar=m=1$.   Equation (\ref{eq:QMeigen}) defines an eigenvalue problem where $\psi(x)$ and $E$ denote the {eigenfunction $f(x,\lambda)$ and eigenvalue $\lambda$ pair}. A  boundary value eigenvalue problem is defined by considering  a  potential function $V(x)$ and boundary conditions of $\psi(x)$.
We assess the performance of the proposed network architecture by solving Eq. (\ref{eq:QMeigen}) for the potential functions of the infinite square well  and the harmonic oscillator, both of which have known analytical solutions.

\subsection{Infinite Square Well}

The infinite square well  problem is characterized by the following potential function:
\begin{align}
\label{eq:infpot}
    V(x) = \begin{cases}
          0 & 0\leq x\leq \ell \\
          \infty & \text{otherwise}
       \end{cases},
\end{align}
 where  the length of the well is set to $\ell=1$. The exact  eigenfunctions and eigenvalues read
\begin{align}
\label{eq:infefunc}
\psi_n(x) = \begin{cases}
          \sqrt{{2}} \sin({n\pi}x) & 0\leq x\leq 1 \\
          0 & \text{otherwise}
       \end{cases},
\quad \quad       E_n = \frac{n^2\pi^2}{2},
\end{align}
where $n$ is a positive integer and  indicates different solutions.
The eigenfunctions are strictly zero outside of the well, implying the boundary conditions $\psi(0)=\psi(1)=0$.  The Eqs. (\ref{eq:parSol}) and (\ref{eq:parametric}) ensure the boundary conditions by setting $\xL=0$, $\xR=1$, and $f_b=0$. The proposed  scanning model is capable of solving for an arbitrary number of the first $n$ states. In Fig. \ref{fig:infiniteplots} we show results up to  $n=3$.
The left panel presents the loss functions of Eqs. (\ref{eq:Loss}) and (\ref{eq:nontriv}) (upper), and the predicted $E$ (lower) during the network optimization.
The scanning algorithm pushes the predicted eigenvalue upwards. The loss falls precipitously when an eigenfunction  is found, and the energy shows plateaus at these exact eigenvalues (indicated by dashed black line) of Eq. (\ref{eq:infefunc}). The loss function in Fig. \ref{fig:infiniteplots} depicts three dips, which correspond to three plateaus in the energy. This behavior gives a physical meaning to the loss function, since by inspecting  $L$ during the training we can draw the eigenstates. The right  panel  shows the extracted $\psi(x)$ (blue) and $E$ (dashed black) at each plateau. Comparing with the exact solutions of Eq. (\ref{eq:infefunc}), the order of magnitude of errors  are $10^{-3}$ and $10^{-4}$  for $\psi_n$ and $E_n$, respectively.
\begin{figure}[h]
    \centering
    \includegraphics[scale=.25]{/inffinal.pdf}
    \caption{Infinite square well: Left panel shows the loss functions and the predicted energy during the training; dashed lines indicate the exact energy levels.  Right plot outlines the predicted eigenfunctions (blue) and eigenvalues (dashed black).  The errors are of the order $10^{-3}$  and $10^{-4}$ for $\psi$ and $E$.}
    \label{fig:infiniteplots}
\end{figure}

\begin{comment}
\subsection{Finite Square Well}

We consider a finite square well potential given by
\begin{align}
\label{eq:finpot}
V(x) = \begin{cases}
          V_0 & x\leq -\frac{L}{2} \\
          0 & -\frac{L}{2}\leq x\leq \frac{L}{2} \\
          V_0 & \frac{L}{2}\leq x
       \end{cases}
\end{align}
and has the following analytical eigenfunctions
\begin{align}
\label{eq:finefunc}
\psi(x) = \begin{cases}
          A\sin({\frac{\sqrt{2mE_n}}{\hbar}}) & -\frac{L}{2}\leq x\leq \frac{L}{2},\; even\;n \\
          A\cos({\frac{\sqrt{2mE_n}}{\hbar}}) & -\frac{L}{2}\leq x\leq \frac{L}{2},\; odd\;n \\
          \pm Be^{\pm \frac{\sqrt{2m(V_0-E_n)}}{\hbar}x} & otherwise
       \end{cases}
\end{align}
Although similar to the infinite square well, it differs in that the eigenfunctions can have non-zero values in the regions with finite potential, a phenomenon called quantum tunneling. Therefore, unlike the infinite case, the eigenfunction solutions should not be fixed to 0 at the boundary walls. However, solutions with eigenvalues smaller than the potential cannot exist in a stable form outside the well, so we  expect the eigenfunction to approach 0 in regions far from the well. This boundary condition at infinity is approximated by parameterizing our eigenfunction solutions to be 0 at the boundaries $x = \pm 6L$.

\begin{figure}[h]
    \centering
    \includegraphics[width=\linewidth]{Latest/finplots_final.png}
    \caption{Training results from the finite square well experiment. Top left shows the various loss terms . Bottom left plots the history of the model's eigenvalue. Middle plot shows the eigenfunctions found by the model. The three plots on the right are the errors of each eigensolution to the known true solution.}
    \label{fig:finiteplots}
\end{figure}

Figure \ref{fig:finiteplots} reveals that our model is able to capturing the quantum tunneling effect. From the middle plot, where the vertical cyan lines depict the well's boundaries, we can see that the non-zero (but exponentially decreasing) eigenfunction is found outside the well. As with the infinite square well, the loss plot depicts three "dips", which correspond to three plateaus in the eigenvalue history of the model during training. These are the points during scanning that our model finds an eigensolution. The error plots again demonstrate the accuracy of our method.
\end{comment}

\subsection{Quantum Harmonic Oscillator}

The harmonic oscillator is characterized by the quadratic potential function
\begin{align}
\label{eq:harmpot}
V(x) = \frac{1}{2}kx^2,
\end{align}
where $k$ is the force constant and is considered to be $k=4$. The exact solutions for the eigenfunctions and energies are  given in terms of Hermite polynomials $H_n$ as
\begin{align}
\label{eq:harmsols}
\psi_n(x) = \frac{1}{\sqrt{2^n n!}}  \frac{e^{-\frac{ x^2}{2}}}{\pi^{1/4}}     H_n\left( x \right),
\quad \quad   E_n =   n + \frac{1}{2} .
\end{align}

The boundary conditions for the quantum oscillator problem  dictate the wavefunction to vanish at infinity, that is, $\psi(-\infty)=\psi(\infty)=0$. In numerical methods,  infinity is assumed to be a large number compared to the potential dimensions. We adopt the same approach and consider the boundary conditions $\psi(-6)=\psi(6)=0$. Thus, Eqs.  (\ref{eq:parSol}) and (\ref{eq:parametric}) ensure the boundary conditions by setting $\xL=-6$, $\xR=6$, and $f_b=0$.

\begin{figure}[h]
    \centering
    \includegraphics[scale=.26]{/harfinal.pdf}
    \caption{Quantum harmonic oscillator: Top left shows the various loss terms. Bottom left plots the history of the predicted energy. Right plot shows the eigenfunctions found by the model; the red dotted line outlines the potential. The order of magnitude of errors are $10^{-2}$ for $\psi$ and $10^{-2}$ for $E$.}
    \label{fig:harmonicplots}
\end{figure}

The proposed scanning neural network method is employed to discover the first three eigenstates for the quantum harmonic oscillator. The left panel in Fig. \ref{fig:harmonicplots} shows the drops in the total loss (blue line) that correspond to plateaus in the eigenvalue during training (green), indicating that an eigenvalue has been found. The predicted $\psi_n(x)$ and $E_n$ are presented in the right panel in Fig.  \ref{fig:harmonicplots} by solid blue and dashed black lines, respectively, while the red dotted curve outlines the potential energy.

\section{Conclusion}
In recent years, there has been a growing interest in the application of  neural networks to study differential equations. In this work, we introduced a neural network  that is capable  of discovering  eigenvalues and eigenfunctions for boundary conditioned differential eigenvalue problems. The obtained solutions identically satisfy the given boundary conditions. A scanning mechanism allows the network to find an arbitrary number of eigenvalues and associated eigenfunctions.
Inspecting the loss function during training allows one to draw the eigenstates,  providing  a physical meaning to the loss function.
The optimization solely depends on the network's predictions, consisting of an unsupervised learning method.  We  demonstrated the capability of the proposed architecture by solving  the   infinite well and harmonic oscillator quantum problems.

\begin{ack}
The authors would like to acknowledge fruitful discussions with Dr. David Sondak.
\end{ack}

\bibliographystyle{plain}


\section{Introduction}
\label{sec:intro}

A knot $K\subset S^3$ is \emph{alternating} if it admits a diagram
with the property that when traversing through the diagram, we alternate
between over- and under-crossings. (An intrinsic definition of alternating
knots have been recently found by Greene and Howie~\cite{Josh, Howie}.)
A diagram of $K$ partitions the plane into domains (the connected components
of the complemet of the projection), and the alternating property
can be rephrased by saying that on the boundary of each domain each edge
connects an under-crossing with an over-crossing. Indeed, this observation
provides a way to measure how far a knot is from being alternating.
We introduce the following definition:

\begin{defn}
  Suppose that $D$ is the diagram of a given knot $K\subset S^3$.
  A domain $d$ of $D$ is \defin{good} if any edge
  on the boundary of $d$
  connects an over- and an undercrossing. The domain $d$ is
  \defin{bad} if it is not good. The number of bad domains of the diagram
  $D$ is denoted by $B(D)$.
\end{defn}

Obviously the diagram $D$ is alternating iff $B(D)=0$. Indeed,
by taking
\[
\beta (K)=\min \{ B(D)\mid D\ \mbox{diagram\ for\ }\ K\}
\]
we get a knot invariant, which satisfies that $\beta (K)=0$ if and
only if $K$ is an alternating knot. As it is typical for knot
invariants given by minima of quantites over all diagrams, it is easy
to find an upper bound on $\beta (K)$ (by determining $B(D)$ for a
diagram of $K$), but it is harder to actually compute its value.

As it turns out, knot Floer homology provides a lower bound for $\beta
(K)$ through the \emph{thickness} of $K$. Recall that ${\widehat {\rm {HFK}}} (K)$, the
hat-version of knot Floer homology of $K$ is a finite dimensional
bigraded vector space over the field ${\mathbb {F}}$ of two elements. By
collapsing the Maslov and Alexander gradings $M$ and $A$ on ${\widehat {\rm {HFK}}}
(K)$ to $\delta =A-M$ we get a graded vector space ${\widehat {\rm {HFK}}}
^{\delta}(K)$. The thickness $th(K)$ of $K$ is the largest possible
difference of $\delta$-gradings of two homogeneous (nonzero) elements
of this vector space. It is known that for an alternating knot $K$ the
$\delta$-graded Floer homology is in a single $\delta$-grading
(determined by the signature of the knot), hence if $K$ is
alternating, then $th(K)=0$. (Knots satisfying $th(K)=0$ are called
\emph{thin} knots, hence alternating knots are thin.)

With this definition in place, the main result of this paper is
\begin{thm}\label{thm:estimate}
Suppose that $K\subset S^3$ is a non-alternating knot. Then
    \begin{equation}\label{eq:ThicknessBadness}
    th(K)\leq \frac{1}{2}\beta (K)-1.
    \end{equation}
\end{thm}

While the thickness of $K$ can be used to estimate how non-alternating
$K$ is, the formula of Equation~\eqref{eq:ThicknessBadness} also can
be used to estimate $th(K)$ by finding appropriate diagrams of $K$.
In particular, the formula can be applied to show

\begin{cor}(Lowrance~\cite{Lowrance})
  \label{cor:Monte}
  Suppose that $K$ is a Montesinos knot. Then, $th(K)\leq 1$.
\end{cor}

\begin{rem}
  A quantity similar to $\beta (K)$ has been introduced by Turaev~\cite{Turaev},
  now called the \defin{Turaev genus} $g_T(K)$. An inequality similar to
  Inequality~\eqref{eq:ThicknessBadness} for the Turaev genus and the
  (knot Floer) thickness $th(K)$ was shown by Lowrence in~\cite{Lowrance}.
  As the Turaev genus of non-alternating Montesinos knots is known
  to be equal to 1~\cite{ChampKof, FiveAuthors},
  our Corollary~\ref{cor:Monte} also follows from \cite{Lowrance}.
\end{rem}

The formula of Inequality~\eqref{eq:ThicknessBadness} can be used in a
further way: by a recent result of Zibrowius~\cite{Zib} mutation does not
change ${\widehat {\rm {HFK}}} ^{\delta}(K)$, hence leaves $th(K)$ unchanged. Consequently,
besides isotopies we can change a diagram by mutations to get better estimates
for $th(K)$ through $B(D)$ for a diagram $D$ of a mutant.

The paper is organized as follows. In Section~\ref{sec:KnotFloer} we
recall basics of knot Floer homology and prove the theorem stated
above. In Section~\ref{sec:Montesinos} we give the details of the
proof of Corollary~\ref{cor:Monte}, and finally in
Section~\ref{sec:further} we list some further properties and
questions regarding $\beta$.

{\bf {Acknowledgements}}: AS was partially supported by the
\emph{\'Elvonal (Frontier) project} of the NKFIH (KKP126683).  ZSz was
partially supported by NSF Grants DMS-1606571 and DMS-1904628. We
would like to thank Jen Hom and Tye Lidman for a motivating
discussion.

\section{The knot Floer homology thickness of knots}
\label{sec:KnotFloer}

Suppose that $V=\sum _aV_a$ is a finite dimensional graded vector space,
where $V_a\subset V$ is the subspace of homogeneous elements of
grading $a\in \mathbb R$. The \emph{thickness} $th(V)$ of $V$ is by definition the
largest possible difference between gradings of (non-zero) homogeneous elements:
\[
th(V)=\max \{ a\in \mathbb R \mid V_a\neq 0\}-\min \{a \in \mathbb R  \mid V_a\neq 0\}.
\]

Suppose now that the graded vector space $V$ is endowed with a
boundary operator $\partial$ of degree $1$; then the homology
$H(V,\partial )$ also admits a natural grading from the grading of $V$.
As $H(V, \partial )$ is the quotient of a subspace of $V$,
it is easy to see that
\[
th(H(V,\partial ))\leq th (V).
\]

The hat version of knot Floer homology (over the field ${\mathbb {F}}$ of two
elements) of a knot $K\subset S^3$ is a finite dimensional bigraded
vector space ${\widehat {\rm {HFK}}} (K)=\sum _{M,A} {\widehat {\rm {HFK}}} _M (K, A)$. By collapsing
the two gradings to $\delta =A-M$, we get the $\delta$-graded
invariant ${{\widehat {\rm {HFK}}} }^{\delta}(K)$. The thickness of ${{\widehat {\rm {HFK}}}
}^{\delta}(K)$ is by definition the thickness $th(K)$ of $K$.

Knot Floer homology is defined as the homology of a chain complex
we can associate to a diagram of the knot (and some further choices).
Indeed, for a given  diagram $D$ of a knot $K$ fix a marking, that is a
point of $D$ which is not a crossing.
Consider the bigraded vector space $C_{D,p}$ (graded by the Alexander
and the Maslov gradings $A$ and $M$) associated to the marked
diagram $(D,p)$, which is  generated over ${\mathbb {F}}$ by the
Kauffman states of the marked diagram, a concept which we recall below.

Suppose that for the marked diagram $(D, p)$ of the knot $K$
the set of crossings is denoted by $Cr(D)$, the set of domains by
$Dom(D)$, and $Dom_p(D)$ denotes the set of those domains which do not contain
$p$ on their boundary. A \emph{Kauffman state} $\kappa$ is a bijection
$\kappa \colon Cr (D)\to Dom_p(D)$ with the property that for a
crossing $c\in Cr (D)$ the value $\kappa (c)$ is one of the (at most
four) domains meeting at $c$. The Alexander, Maslov and
$\delta$-gradings of a Kauffman state are computed by summing the
local contributions at each crossing, as given by the diagrams of
Figure~\ref{fig:Contributions}.
\begin{figure}
\centering
\includegraphics{Contributions.eps}
\caption{{\bf The local contributions for $A, M$ and $\delta$ at a
    crossing.}  The Kauffman state distinguishes a corner at the
  crossing, and we take the value in that corner as a contribution of
  the crossing to $A,M$ or $\delta$ of the Kauffman state at hand.}
\label{fig:Contributions}
\end{figure}

According to \cite{OSzAlternating} there is a boundary map $\partial
\colon C_{D,p}\to C_{D,p}$
of bidegree $(-1,0)$ (in the bigrading $(M,A)$)
with the property that $H(C_{D,p}, \partial
)={\widehat {\rm {HFK}}} (K)$ is isomorphic to the knot Floer homology of $K$
(as a bigraded vector space). By collapsing the two
gradings $A$ and $M$ to $\delta =A-M$, we get the graded vector spaces
$(C_{D,p}^{\delta},\partial )$ and its homology ${{\widehat {\rm {HFK}}}
}^{\delta}(K)$.  As ${\widehat {\rm {HFK}}}^{\delta }(K)$ is the quotient of a subspace
of $C^{\delta}_{D,P}$, we have that $th({{\widehat {\rm {HFK}}} }^{\delta}(K)) \leq
th(C_{D,p}^{\delta}, \partial)$.

\begin{prop}
  \label{prop:ThicknessBadDomains}
  Suppose that $D$ is a diagram of the knot $K$.  If $D$ is not an
  alternating diagram, then
  \[
 th(C^{\delta}_{D,p})\leq \frac{1}{2}B(D)-1.
\]
  \end{prop}
\begin{proof}
Fix a marked point $p$ on $D$, and
  consider the $\delta$-graded chain complex
  $(C^{\delta}_{D,p},\partial )$ generated by the Kauffman states of $(D,p)$.

  The $\delta$-grading at a positive crossing is either 0 or $-\frac{1}{2}$,
and at a negative crossing it is either 0 or $\frac{1}{2}$.
So we can express the $\delta$-grading of a
  Kauffman state $\kappa$ as the sum
  \[
  -\frac{1}{4}{\rm {wr}}(D)+\sum _{c\in Cr}f(\kappa (c)),
  \]
  where ${\rm {wr}}$ is the writhe of the diagram, and $f$ is a function
  on the Kauffman corners, which is either $\frac{1}{4}$ or $-\frac{1}{4}$
  (depending on the chosen quadrant at the crossing $c$).

  Simple computation shows that for a good domain each corner in the
  domain gives the same $f$-value, hence for different Kauffman states
  the contributions from this particular domain are the same. This is
  no longer true for a bad domain, but the difference of two
  contributions is at most $\frac{1}{2}$.  When determining the possible
  maximum of $\delta (x)-\delta (x')$ for two homogeneous elements
  $x,x'\in C_{D,p}^{\delta}$, the contributions from the writhe cancel,
  and so do the contributions from good domains, while bad domains
  contribute at most $\frac{1}{2}$. This shows that
  $th(C_{D,p}^{\delta})\leq \frac{1}{2}B(D)$.

  By assumption $D$ is not alternating, hence there is a bad domain,
  with an edge showing that it is bad. Choose the marking $p$ on such an
  edge.  Since this edge guarantees that the two domains having it
 on their boundary are both bad, while these two bad domains do not
  get Kauffman corners, we get that $th(C_{D,p})$ is bounded by
  $\frac{1}{2}(B(D)-2)=\frac{1}{2}B(D)-1$, concluding the proof.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{thm:estimate}]
  Suppose that $K$ is not alternating. Then any diagram $D$ of
  $K$ is non-alternating, hence we have that
\[
  th(K)\leq th(C_{D,p}^{\delta})\leq \frac{1}{2}B(D)-1.
  \]
  Since $\beta (K)$ is computed from the minimum of the right-hand
  side of this inequality, the proof follows at once.
\end{proof}

\section{Montesinos knots}
\label{sec:Montesinos}

Montesinos knots are straighforward generalizations of pretzel knots;
a diagram involving rational tangles defining the Montesinos
knot $M(r_1, \ldots , r_n)$ is shown by Figure~\ref{fig:monte}.
(A box with a rational number $r_i$ in it symbolizes the
tangle shown by Figure~\ref{fig:tangle}.)
\begin{figure}
\centering
\includegraphics{monte.eps}
\caption{{\bf The Montesinos knot $M(r_1, \ldots , r_n)$.}  The box
  containing $r_i$ denotes the algebraic tangle determined by the
  rational number $r_i=\frac{\beta _i}{\alpha _i}$
  (cf. Figure~\ref{fig:tangle}).  In order to have a knot, at most one
  of $\alpha _i$ can be even.}
\label{fig:monte}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=5cm]{tangle.eps}
\caption{{\bf The rational tangle corresponding to $r\in \mathbb Q$.}  Here
  the boxes with $c_i\in \mathbb Z$ on the right denote $\vert c_i\vert $
  half twists (right handed for positive, left handed for negative
  $c_i$). Depending on the parity of $n$ (the numbre of $c_i$) we have
  two different finishing forms. The rational number $r$ determines
  the coefficients $c_i$ through its continued fraction expansion. The
  tangle is alternating (as part of a knot or link) if $c_i$ alternate
  in signs.}
\label{fig:tangle}
\end{figure}
We allow any of the $r_i$ to be equal to $\pm 1$. Notice that the
order of $(r_1, \ldots , r_n)$ is important; those $r_i$ which are equal
to $\pm 1$ can be commuted with any other parameter through a simple
isotopy of the diagram.

\begin{lem}\label{lem:isotopy}
  Consider the diagram of the Montesinos knot $M(r_1, \ldots , r_n)$ given
  by Figure~\ref{fig:monte}.
  It can be isotoped to a diagram with at most four bad domains.
\end{lem}

\begin{proof}
  Recall that a rational tangle has the form given by Figure~\ref{fig:tangle}.
Adapting the isotopies described in \cite{Goldman-Kauffman}, we can
achieve that all tangles are alternating, hence the potentially bad
domains are the ones between the tangles, together with the central
and the unbounded domains.  The number of bad domains
between the tangles can be reduced by the following observation.  The
domain between two tangles is bad if the first coefficients $c_1^1$
and $c_1^2$ of the two rational numbers determining the tangles have
opposite signs, say $c_1^1>0$ and $c_1^2<0$. Then by Reidemestier-II
moves we can introduce cancelling twistings, as shown by
Figure~\ref{fig:cancel}, and then commute the first twisting (in the figure
given by the box with $2$ in it) between the first and second tangles
of the Montesinos knot.
\begin{figure}
\centering
\includegraphics[width=5cm]{cancel.eps}
\caption{\bf The introduction of cancelling twists to turn domains
    between tangles good.}
\label{fig:cancel}
\end{figure}
All domains between the boxes will become good, except the ones
connecting the first tangle with the newly introduced twists and the
second tangle also connecting it with the newly introduced twists.
After these alterations make sure that (by the adaptation of
\cite{Goldman-Kauffman}) all tangles are isotoped to be alternating.
In total the new diagram then has four bad domains, concluding the
proof.
\end{proof}

\begin{proof}[Proof of Corollary~\ref{cor:Monte}]
  For a Montesinos knot $M(r_1, \ldots , r_n)$ an approrpiate isotopy
  of the diagram of Figure~\ref{fig:monte} (as given by
  Lemma~\ref{lem:isotopy}) gives a diagram with at most four bad
  domains. The application of Theorem~\ref{thm:estimate} concludes the
  argument.
\end{proof}

\begin{rem}
  Using the mutation invariance of $th(K)$, Lemma~\ref{lem:isotopy}
  can be avoided: by mutations any Montesinos knot $M(r_1, \ldots ,
  r_n)$ can be moved to $M(q_1, \ldots , q_n)$ with the same rational
  parameters in a different order so that $q_i$ and $q_{i+1}$ have the
  same sign with at most one exception. Isotoping the diagram so that
  the tangles are alternating, the mutated diagram then has at most 4
  diagrams. Using the result of \cite{Zib} then the corollary follows
  as before.
  \end{rem}

\section{Further properties}
\label{sec:further}
It is a standard fact that the knot Floer homology of the connected
sum of two knots is the tensor product of the knot Floer homologies:
\[
{\widehat {\rm {HFK}}} (K_1\# K_2)\cong {\widehat {\rm {HFK}}} (K_1)\otimes {\widehat {\rm {HFK}}} (K_2).
\]
From this (bigraded) isomorphism it follows that
\[
th(K_1\# K_2)=th(K_1)+th(K_2).
\]
The behaviour of $\beta (K)$ is less clear under connected summing.
Suppose that $K_1, K_2$ are both non-alternating knots. By taking the
connected sum of two diagrams $D_1,D_2$ for these knots at bad edges
(i.e. arcs on the boundary of bad domains verifying that the domains
are bad) we get that
\[
B(D_1\# D_2)=B(D_1)+B(D_2)-2,
\]
immediately implying that
\[
\beta (K_1\# K_2)\leq \beta (K_1)+\beta (K_2)-2.
\]
Motivated by the equality for the thickness $th$, we conjecture
\begin{conj}
  If $K_1, K_2$ are two non-alternating knots, then
  \[
  \beta (K_1\# K_2)= \beta (K_1)+\beta (K_2)-2.
  \]
\end{conj}

\bigskip

It is not hard to find knot diagrams for which
Inequality~\eqref{eq:ThicknessBadness} is sharp. Indeed, the standard
diagram of the pretzel knot $P(-3,5,5)$ admits four bad domains (see
Figure~\ref{fig:equality}(a)), while an explicite calculation of
${\widehat {\rm {HFK}}} (P(-3,5,5))$ shows that $th(P(-3,5,5))=1$. Consider the
$n$-fold connected sum $K_n=\#_n P(-3,5,5)$;  connect summing the
diagrams at bad edges (in the above sense) we get a sequence of knots
$K_n$ and diagrams $D_n$ for them with the properties that $th
(K_n)=n$ and $B(D_n)=2n+2$, cf. Figure~\ref{fig:equality}(b).
The non-alternating knots $K_n$ then satisfy
$n=th(K_n)=\frac{1}{2}\beta (K_n)-1$.
\begin{figure}
\centering
\includegraphics[width=5cm]{equality.eps}
\caption{ In (a) he pretzel knot $P(-3,,5,5)$ is shown. The $B$ symbols
  signify the bad domains. (A box containing
  the integer $n$ denotes $\vert n \vert $ half twists, right handed for $n>0$
  and left handed for $n<0$.) In (b) we provide a diagram of $K_n$, where the
connected sum is taken at bad domains.}
\label{fig:equality}
\end{figure}


\section{Introduction}

Memristors have been proposed as new electronic devices with promising analog/digital and neuromorphic applications
\cite{shin2010memristor,zidan2014compensated,fouda2019spiking,wu2019enhanced,karimi2019novel,fouda2014memristor}. Their characteristic behavior is a pinched hysteresis loop in the current-voltage plane \cite{ascoli2015nonlinear}. Many circuits (known as emulators) have been introduced to mimic the memristor's behavior  \cite{ranjan2018high,petrovic2019tunable,elwakil2013simple,alharbi2017electrical,fouda2014simple}. Although this behavior has been shown to exist in other nonlinear devices such as a non-linear inductor or a non-linear capacitor \cite{fouda2015pinched,corinto2012memristive,sadecki2019analysis}, and that its appearance is linked to satisfying the necessary conditions of the theory of Lissajous figures  \cite{maundy2019correlation, maundy2018simple, weilnhammer2016perceptual}, it is still mistakenly attributed  to memristors only. There are actually doubts about the uniqueness of the memristor as a fundamental device which have been raised by several researchers (see \cite{abraham2018case} and the references therein) but there is no doubt that it is a dynamic and nonlinear device \cite{riaza2014comment,elwakil2013simple, non}. Without being nonlinear, the frequency-doubling mechanism mandated by the theory of Lissajous figures, which is essential to create a pinched loop (as explained in detail in \cite{maundy2019correlation} and most recently in \cite{MAJZOUB}) cannot be obtained. It is important to note that even if the memristor was not a fundamental device, it may still be a useful nonlinear device. Therefore, semi-conductor devices labeled as memristors continue to be fabricated in different technologies and materials \cite{babacan2018fabrication, bio} and are now commercially available from some vendors (see for example www.knowm.com).

Any memristive system is generally defined as
\cite{chua1976memristive,radwan2015mathematical}

\begin{align}
y=g(x,u,t)u \nonumber\\
\dot{x}=f(x,u,t) \nonumber
\end{align}
where $u$ and $y$ are the input and output of the system and $x$ is the state variable. To identify a device as a memristor, it has to have three fingerprints
\cite{adhikari2013three,radwan2015mathematical}
: 1) a pinched loop  in current-voltage plane with
2)  a loop area that decreases monotonically with increasing frequency and 3) shrink to a single-valued function when the frequency tends to infinity. This  does not mean that pinched hysteresis loops cannot be observed in other devices or circuits (with increasing or decreasing loop area).

In this paper, we provide  two examples of  simple fundamental circuits, namely  nonlinear series and parallel resonators, that can exhibit the pinched hysteresis loop behavior. Mathematical models are derived and numerically validated with a diode-type asymmetric switching nonlinearity and also with an odd-symmetric nonlinearity formed of two anti-parallel diodes. Experimental results are  provided and confirm the theory. The formation of a pinched loop from these circuits confirms our previous assertions \cite{fouda2015pinched,maundy2019correlation} that:
\begin{enumerate}
\item pinched hystersis is not a finger-print of memristors only,
\item  pinched loops with lobe area widening (rather than declining) with increased frequency are possible and
\item the existence of a nonlinearity is essential for generating this behavior. This is because the theory of Lissajous figures necessitates the generation of a second harmonic frequency in the electrical current $i(t)$ when the applied voltage on the device $v(t)$ has only one fundamental frequency (see \cite{MAJZOUB} for a more detailed explanation).
\end{enumerate}

We further demonstrate the use of the generated pinched loop behavior in realizing digital logic AND/OR gates. With this, we aim to show that research in the possible applications of pinched loops \cite{app1, app2} can be conducted without restricting this behavior to memristors.

\section{Nonlinear resonators}

\subsection{Series RLC Resonator}

A simple nonlinear resonator circuit is shown in Fig. \ref{fig:sercit} composed of a series RLC resonance circuit interrupted by a voltage-controlled switching nonlinearity (denoted $f(x)$) placed across the capacitor. The circuit is excited by a sinusoidal signal $v_{in}(t)=V_{DC}+V_A\sin(\omega t)$. In its simplest form, the switching device can be a single diode and in this case the circuit is described by the equation set

\begin{subequations}
\begin{align}
C\frac{dv_{C}}{dt}=i_{L}-i_{f},\\
L\frac{di_{L}}{dt}=v_{in} -Ri_{L}-v_{C}
\end{align}
\end{subequations}
where $v_C$ and $i_L$ are respectively the voltage and current across the capacitor and in the inductor while $i_f$ is the nonlinear device current given by

\begin{figure}[h]
\centering{\includegraphics[width=60mm]{Schematic.pdf}}
\caption{Nonlinear series resonator circuit with diode-based nonlinear switches.}
\label{fig:sercit}
\end{figure}

\begin{equation}
i_f=i_{D}=\left\{ \begin{array}{cc}
(v_{C}-v_{\gamma})/R_{D} & \:\:v_{C}\geq v_{\gamma}\\
0 & \:\:v_{C}<v_{\gamma}
\end{array}\right.
\end{equation}
where $R_D$ is the diode switch-on resistance while $v_{\gamma}$ is the switch-on voltage (approximately $0.6 V$ for silicon devices). Introducing the dimensionless variables: $x=v_{C}/v_{\gamma},y=Ri_{L}/v_{\gamma},a=R/R_{D},b=R^2C/L, A=V_{DC}/v_{\gamma}, B=V_A/v_{\gamma}$, a normalized frequency $\omega_n=\omega RC$ and a normalized time $t_n=t/RC$, the above equations transform into
\begin{subequations}
\begin{align}
\frac{dx}{dt_n}=y-f(x),\\
\frac{dy}{dt_n}=b\left(g(t_n )-y-x\right)
\end{align}
\end{subequations}
with $g(t_n)=A+B\sin(\omega_n t_n)$  being the sinusoidal excitation and $f(x)$ being the nonlinearity. In the case of the switching-diode, $f(x)$ is asymmetric and  given by

\begin{equation}
f(x)=\left\{ \begin{array}{cc}
a(x-1) & \;x\geq1\\
0 & \;x<1
\end{array}\right.
\end{equation}
Numerical simulations of this model were performed using the ode45 solver in MATLAB with $A=1,b=0.2,B=2$ for different values of $a$ and $\omega_n$ as shown in Figs. \ref{fig:serNum}(a) and 2(b) which plot the observed pinched loop in the $g(t_n)-y(t_n)$ plane; i.e.  in the $v_i-i_L$ plane of the circuit. Therefore, this pinched-loop represents an impedance analogous to the "memristance". It is noted from these simulations that the lobe area increases with increasing the value of the parameter $a$  and also with increasing the frequency $\omega_n$. Note that $\omega_n$ can also be written as $\omega_n=\omega\cdot BW/\omega_{r}^{2}$ where $\omega_r=1/\sqrt{LC}$ is the resonance frequency and $BW=R/L$ is the resonance bandwidth. Hence, satisfying a certain value of $\omega_n$ is achieved not only by changing the applied excitation signal frequency $\omega$ but also by adjusting the resonance bandwidth and resonance frequency.

\begin{figure}[h]
 \centering{\includegraphics[width=130mm]{Different.pdf}}
 \caption{Numerical simulation results of the nonlinear model of the series resonator circuit with (a), (b) a single diode-nonlinearity  and (c), (d) with two anti-parallel diodes.  The effect of  different values of parameter $a$ at $\omega_n=0.1$ are shown in  (a) and (c) and the effect of different values of $\omega_n$ at $a=1$ are shown in  (b) and (d)}
\label{fig:serNum}
\end{figure}

Based on the detailed analysis and explanation of \cite{maundy2019correlation}, the appearance of the pinched loop in this nonlinear resonator is a result of satisfying the conditions of the theory of Lissajous figures. In particular, the nonlinearity $f(x)$ is responsible for generating a strong second-harmonic  in the current signal $y(t)$ with a normalized frequency $2\omega_n$. The existence of this harmonic can easily be verified from the  FFT of $y(t)$. So long as the phase shift between this harmonic and the fundamental frequency component remains less than $\pi/2$,  the pinched-loop will be observable \cite{maundy2019correlation}. The location of the pinch point can be found  using the general equations derived in \cite{maundy2019correlation}.  When the diode-type nonlinearity is replaced with an odd-symmetric one, such as a cubic-type nonlinearity  or the nonlinearity resulting from two anti-parallel diodes given by (5) (see appendix), two pinched-points are obtained in this case as shown in Figs. \ref{fig:serNum}(c) and (d)  for the parameter values $A=0,b=0.2$  and $B=2$.A pure sinusoid  is used in this case for excitation ($A=0$).

\begin{equation}
f(x)=\left\{ \begin{array}{cc}
a(x-1) & \;x\geq1\\
0 & \;-1<x<1\\
a(x+1) & \;x\leq -1
\end{array}\right.
\end{equation}

Shaping of the pinched-loop is thus clearly related to the choice of the  function $f(x)$. In particular, the odd-symmetric nonlinearity generates a third-order harmonic from the applied fundamental frequency of the excitation voltage. The total number of pinch points in a Lissajous figure formed from two sinusoids with frequencies $p \cdot \omega$ and $q \cdot \omega$ is known to equal $q(p-1)+p(q-1)$ which is equal to 2 when $p=1$ and $q=3$.

\subsection{ Parallel RLC Resonator}

Another  nonlinear resonator circuit is shown in Fig. \ref{fig:parcit} composed of a parallel RLC resonance circuit interrupted by nonlinearity placed in series  with the inductor. $C_p$ is a parasitic capacitor that must be included in the analysis since the non-linearity is voltage-controlled. In this case the circuit is described by the equation set

\begin{figure}[h]
\centering{\includegraphics[width=45mm]{par.JPG}}
\caption{Nonlinear parallel resonator circuit}
\label{fig:parcit}
\end{figure}

\begin{subequations}
\begin{align}
RC\frac{dv_{C}}{dt}=v_{in}-v_C-i_L R, \\
L \frac{di_L}{dt}=v_c-v_{c_p},\\
C_p \frac{v_{c_p}}{dt}=i_L-i_f
 \end{align}
\end{subequations}

With the same dimensionless variables introduced in the previous subsection, the above equations can be written as

\begin{subequations}
\begin{align}
\frac{dx}{dt_n}=g(t_n )-x-y, \\
\frac{dy}{dt_n}={b}(x-z)\\
\epsilon \frac{dz}{dt_n}=y-f(z)
\end{align}
\end{subequations}
where $z=i_D R/v_\gamma$, $\epsilon=C_p/C$ and $f(z)$ is the switching nonlinearity function given by (4) or (5). When $\epsilon$ tends to zero, (7c) yields $y=f(z)$.

\begin{figure}[h]
\centering{\includegraphics[width=0.8\linewidth]{Different_parallel}}
\caption{Numerical simulation results of the nonlinear model of the parallel resonator circuit with a single diode-nonlinearity and with two anti-parallel diodes: (a), (c) at $\omega_n=0.1$ for different values of $a$ and (b), (d) {at $a=1$ } for different values of $\omega_n$}
\label{fig:parNum}
\end{figure}

Figure \ref{fig:parNum} shows the numerical simulation of (7) with $A=1,b=0.2,B=2$  and $\epsilon=0.01$. The area enclosed inside the lobes increases with increasing the frequency $\omega_n$ but decreases with increasing the value of parameter $a$.

\section{Experimental results}
The nonlinear series resonator circuit was experimentally tested with discrete components $L=1mH$, $C=5.6 nF$, $R=10k \Omega$ and using a 1N914 diode. Figure \ref{fig:serexp} shows the observed pinched loop when the excitation voltage has an offset voltage $V_{DC}=1V$, an amplitude of $V_A=2V$ and at different frequencies. The current in the resistor was measured using a differential probe.\\
At $100Hz$ (see Fig. \ref{fig:serexp}(a)) it is seen that the $v_i-i_L$ portrait resembles the diode characteristics as expected. As frequency increases, the pinched loop formation is observed in Figs. \ref{fig:serexp}(b), (c) and (d). The voltage and current  waveforms corresponding to the portrait of Fig. \ref{fig:serexp}(d) are shown in Fig. \ref{fig:serexptime}. Adding another diode with opposite polarity in parallel to the diode in the circuit, a loop with two pinch points is seen in Fig. \ref{fig:serexpfft}(a) as predicted numerically. The applied voltage in this case  had no offset value ($V_{DC}=0V$), an amplitude of $V_A=7V$ and a $1kHz$ frequency. The corresponding power spectrum of the current signal $i_L$ is shown in Fig. \ref{fig:serexpfft}(b) from which it is visible that strong odd harmonics have been generated. However, the two pinched points are created because the phase angle between the fundamental  signal at $1kHz$ and the generated third-order harmonic at $3kHz$ is less than $\pi/2$ satisfying the Lissajous conditions \cite{maundy2019correlation}. The remaining harmonics contribute to the shape of the loop but not to the creation of the pinch-points.

\begin{figure}[h]
\centering{\includegraphics[width=70mm]{serexp}}
\caption{Results from an experimental setup: (a)-(d) formation of the pinched loop at different frequencies (Xaxis $v_{i}(t)$: 0.6V/div, Yaxis: $i_{L}(t)$: 0.3V/div).}
\label{fig:serexp}
\end{figure}

\begin{figure}[h]
\centering{\includegraphics[width=60mm]{serexptime}}
\caption{Time domain waveforms at $1.5kHz$}
\label{fig:serexptime}
\end{figure}

\begin{figure}[h]
\centering{\includegraphics[width=85mm]{2-diodes.pdf}}
\caption{Experimental results using anti-parallel diodes: (a) loop with two pinch points  (Xaxis $v_{i}(t)$: 2V/div, Yaxis: $i_{L}(t)$: 0.65V/div). (b) power spectrum of $i_L(t)$ (Xaxis: 1.25kHz/div, Yaxis: 10dB/div)}
\label{fig:serexpfft}
\end{figure}

\section{Application and Comparison with a Commercial Memristor Device}
Digital application of AND/OR gates are used to compare their functionality with pinched loops created (i) from the series resonator described above and (ii) from a KNOWM commercial memristor device using their discovery kit \cite{knowm}. The gates are built using two memristor devices, as shown in Fig. \ref{fig:OrAND}, driving a capacitive load (C = 10nF).  Figure \ref{fig:OrANDexp} shows the waveform on the capacitor for both the AND/OR gates after using the experimental memristance data of the pinched loop obtained from the nonlinear resonator circuit and that obtained from  a commercial memristor device \cite{knowm}. It is clear from the figure that similar responses are obtained in both cases. Therefore, from an application point of view, what matters is the existence of the pinched loop rather than its origin. On the other hand, resonators are easily realizable on the micro and nano scales \cite{res1, res2} which paves the way for possibly easier implementations of memristors based on the technique presented here. In terms of the non-volatility of these circuits, a procedure to dis-connect the ground terminal connection should be devised in order to maintain the states of the state variables (inductor current and capacitor voltage) un-changed. It is important to mention, that the non-volatility of memristors is not related to  whether they are fundamental  devices or not.

\begin{figure}[h]
\centering{\includegraphics[width=70mm]{ANDOR}}
\caption{Realization of AND/OR gates using two memristors}
\label{fig:OrAND}
\end{figure}

\begin{figure}[h]
\centering{\includegraphics[width=80mm]{MemristorANDOR}}
\caption{Comparison of the AND/OR gate responses using the pinched loop data from the nonlinear series resonator circuit versus the data from a commercial memristor}
\label{fig:OrANDexp}
\end{figure}

\section{Conclusion}
We have verified numerically and experimentally the existence of a dynamic pinched-loop behavior in a simple nonlinear resonator circuits. The key behind generating such a loop is to produce in the electrical current a second-order harmonic (for a single pinch point case) or third-order harmonic (for a two pinch point case) from the applied fundamental frequency of the excitation voltage. When these harmonics have proper phase and amplitude relations with respect to the fundamental \cite{maundy2019correlation}, pinched loops appear. The circuits presented here rely on the  diode nonlinearity as opposed to MOS transistor quadratic nonlinearity \cite{maundy2018simple, babacan2018fabrication} or multiplier-type nonlinearity \cite{elwakil2013simple, complexity}. The proposed systems represent  second-order and third-order dynamical systems which is not the case of  memrisitve systems as defined in \cite{chua1976memristive}.


\section{Introduction.}
Let $D$ be a quadratic discriminant, meaning that $D$ is a nonsquare integer with $D\equiv 0, 1\pmod{4}$. Set $D = 4d+\sigma$, where $\sigma \in \{0,1\}$, and let $\tau = \frac{\sigma+\sqrt{D}}{2}$. It is easy to check that $\tau^2 \in \mathbb{Z}+\mathbb{Z}\tau$, so that
\begin{align*} \mathbb{Z}[\tau] &= \mathbb{Z}+\mathbb{Z}\tau\\
&= \left\{\frac{u+v\sqrt{D}}{2}: u, v \in \mathbb{Z}, u\equiv vD\!\!\!\pmod{2}\right\}. \end{align*}
In what follows, we write $\mathcal{O}_D$ (for ``order of discriminant $D$'') in place of $\mathbb{Z}[\tau]$.

Our aim with this note is to showcase a simple proof of the following criterion for unique factorization in $\mathcal{O}_D$. We remind the reader that if $R$ is a domain then $\pi\in R$ is \emph{irreducible} if $\pi$ is nonzero and not a unit, and if whenever $\pi=\alpha\beta$ with $\alpha,\beta \in R$, either $\alpha$ or $\beta$ is a unit. The element $\pi\in R$ is \emph{prime} if $\pi$ is nonzero and not a unit, and if whenever $\pi\mid \alpha\beta$ (with $\alpha,\beta \in R$) either $\pi\mid\alpha$ or $\pi\mid \beta$; equivalently, a prime is a nonzero element of $R$ for which the principal ideal $(\pi)$ is a prime ideal of $R$. Prime elements are always irreducible; the converse holds in a UFD (unique factorization domain), but not in general.

\begin{thm}\label{thm:main} Suppose that every rational prime number
\begin{equation}\label{eq:pbound} p \le \begin{cases} \sqrt{|D|/3} &\text{if $D < 0$,} \\
\sqrt{D/5} &\text{if $D > 0$}
\end{cases} \end{equation}
that is irreducible in $\mathcal{O}_D$ is also prime in $\mathcal{O}_D$. Then $\mathcal{O}_D$ is a unique factorization domain.
\end{thm}

\subsection{Examples.}
\begin{enumerate}
\item[(i)] [$D=73$] Since $\sqrt{73/5} = 3.8\dots$, the conditions of Theorem \ref{thm:main} concern only the primes $p=2$ and $p=3$. Neither 2 nor 3 is irreducible, since
    \[ 2 = \frac{9+\sqrt{73}}{2} \cdot \frac{9-\sqrt{73}}{2}, \quad \text{while} \quad 3 = (2\sqrt{73}+17)\cdot (2\sqrt{73}-17). \]
    (It is easy to check that all of the factors listed here are nonunits.)
    We conclude that $\mathcal{O}_D=\mathbb{Z}[\frac{1+\sqrt{-73}}{2}]$ is a UFD.

    The number $73$ is not particularly special.\footnote{See \cite{PS19} for a counterpoint to this claim.} It is widely believed that there are infinitely many $D>0$ for which $\mathcal{O}_D$ is a UFD. In fact, Cohen and Lenstra have precise conjectures predicting, for instance, that $\mathcal{O}_p$ is a UFD for 75.44\dots\% of primes $p\equiv 1\pmod{4}$ (see \cite[\S5.10]{cohen93} and \cite{CL84b,CL84a,tRW03}).
\item[(ii)] [$D=-163$] Since $\sqrt{163/3} = 7.3\dots$, we must check $p=2, 3, 5, 7$. As $\tau=\frac{1+\sqrt{-163}}{2}$ is a root of the monic irreducible polynomial $X^2-X+41$, we have that $\mathbb{Z}[\tau] \cong\mathbb{Z}[x]/(X^2-X+41)$. Hence, for each prime $p$,
    \[ \mathbb{Z}[\tau]/(p)  \cong (\mathbb{Z}[X]/(p))/(X^2-X+41) \cong \mathbb{F}_p[x]/(X^2-X+41). \]
     It is straightforward to check that $X^2-X+41$ is irreducible modulo $p$ for each  of $p=2,3,5,7$. (For the odd primes $p$ in this list, it suffices to observe that the discriminant $-163$ of $X^2-X+41$ is a nonsquare mod $p$.) Therefore, $\mathbb{Z}[\tau]/(p)$ is a field, whence $(p)$ is a prime ideal of $\mathcal{O}_D$ and $p$ is a prime element. So  the criterion of Theorem \ref{thm:main} is again satisfied and  $\mathcal{O}_D$ is a UFD. The number $-163$ \emph{is} special; as shown by Heegner, it is the largest (in absolute value) negative $D$ for which $\mathcal{O}_D$ is a UFD (\cite{heegner52}; see also \cite{cox13}).
 \end{enumerate}

We do not claim that Theorem \ref{thm:main} is new. When $\mathcal{O}_D$ is the full collection of algebraic integers inside $\mathbb{Q}(\sqrt{D})$ (the so-called ``maximal order''), basic algebraic number theory says that $\mathcal{O}_D$ is a Dedekind domain with finite class group. Furthermore, results from the geometry of numbers imply that every ideal class is represented by an ideal with norm  bounded by the quantities appearing on the right of \eqref{eq:pbound} (see \cite[Theorem 13.7.10, p.\ 399]{artin11} for $D < 0$ and \cite[Exercise 17, p.\ 300]{cohen93} for $D > 0$). So Theorem \ref{thm:main} follows easily (in this case).

It seems of some interest --- e.g., for the teaching of basic courses in algebra and number theory --- to give a proof of Theorem \ref{thm:main} requiring as little machinery as possible. Several close relatives of Theorem \ref{thm:main} have been proved in the literature without reference to algebraic number theory; see \cite{fendel85,gyarmati83,lanczi65,
popovici57b,popovici57,ramirez16,snyder07,zaupper83,zaupper90}. However, all of these papers either establish results weaker or less complete than  Theorem \ref{thm:main}, or their proofs depend on auxiliary results from the geometry of numbers or the theory of Diophantine approximation. (For example, the beautifully simple method of Ramirez V.\ in \cite{ramirez16} gives a very satisfactory result when $D<0$, but only a partial result for $D>0$.)  Apart from a few easy lemmas concerning the ``norm'' map (see the Notation section below), our proof of Theorem \ref{thm:main} is self-contained, resting only on the commutative ring theory seen in a first graduate algebra course.

\subsection*{Notation.} We let $K$ be the fraction field of $\mathcal{O}_D$, so that $K = \mathbb{Q}(\sqrt{D})$, and we denote conjugation in $K$ with a bar. The \emph{norm} of $\alpha\in K$, denoted $N(\alpha)$, is defined by $N(\alpha) = \alpha\bar{\alpha}$. We recall that $N(\alpha\beta) = N(\alpha)N(\beta)$ for all $\alpha,\beta \in K$, that the norm sends nonzero elements of $\mathcal{O}_D$ to nonzero integers, and that $\alpha$ is a unit of $\mathcal{O}_D$ if and only if $N(\alpha)=\pm 1$. Readers are invited to prove these  results themselves; alternatively, they may consult, e.g., \cite[Chapter 2]{lehman19}.

\section{Proof of Theorem \ref{thm:main}.}

Our proof makes crucial use of the following lemma, which also features in the arguments of \cite{ gyarmati83,popovici57,ramirez16,snyder07,zaupper83,zaupper90}.

\begin{lem}\label{lem:primenorm} Let $\alpha \in \mathcal{O}_D$. If $N(\alpha)=\pm p$, where $p$ is a rational prime, then $\alpha$ is prime in $\mathcal{O}_D$.
\end{lem}

\begin{proof}

Since $\alpha \bar{\alpha}= \pm p$, there is a canonical surjection $\mathcal{O}_D/(p) \twoheadrightarrow \mathcal{O}_D/(\alpha)$. Since $\bar{\alpha}$ is not a unit, the  corresponding kernel is nontrivial (containing, e.g., $\alpha\bmod{p}$). Thus, $\#\mathcal{O}_D/(\alpha)$ is a proper divisor of $\#\mathcal{O}_D/(p) = p^2$. (The last equality comes from noting that $a+b\tau$, for $0\le a, b < p$, form a complete residue system mod $p$.) Since $\alpha$ is not a unit, $\#\mathcal{O}_D/(\alpha) > 1$. Therefore, $\#\mathcal{O}_D/(\alpha)=p$, and so $\mathcal{O}_D/(\alpha) \cong \mathbb{F}_p$. Hence, $(\alpha)$ is a prime (in fact, maximal) ideal of $\mathcal{O}_D$, so that $\alpha$ is prime in $\mathcal{O}_D$.
\end{proof}

We turn now to the proof of Theorem \ref{thm:main}. A simple induction on $|N(\alpha)|$ shows that every nonzero, nonunit $\alpha \in \mathcal{O}_D$ has a factorization into irreducibles. So it remains only to prove uniqueness. We  reduce this (as in \cite{popovici57, snyder07,zaupper90}) to  the following claim.

\begin{claim} Every prime in $\mathbb{Z}$ factors as a product of primes in $\mathcal{O}_D$.
\end{claim}

To see why this suffices, recall that an element with a factorization into primes necessarily has this as its only factorization into irreducibles (up to order and unit factors). This is clear from the usual proof of unique factorization in  a Euclidean domain or PID (compare with the proof of  Proposition 12.2.14(a) in \cite{artin11}). Since every rational integer larger than $1$ factors as a product of rational primes, our claim implies that all those integers factor uniquely in $\mathcal{O}_D$. But this implies that every $\alpha \in \mathcal{O}_D$, not zero and not a unit, also factors uniquely: If $\alpha$ had two factorizations, we could cook up two factorizations of $|N\alpha| = \pm \alpha\bar{\alpha}$ by concatenating our factorizations of $\alpha$ with a fixed factorization of $\pm \bar{\alpha}$.

\begin{proof}[Proof of the claim]  Assuming the claim to be false, let $p$ be the smallest prime for which it fails. Then
\begin{equation}\label{eq:plower} p >
\begin{cases} \sqrt{|D|/3} &\text{when $D<0$}, \\
\sqrt{D/5} &\text{when $D>0$}.
\end{cases}
\end{equation}
Indeed, suppose otherwise. Since $p$ does not factor as a product of primes, it itself is not prime. But then the hypothesis of Theorem \ref{thm:main} tells us that $p$ factors nontrivially in $\mathcal{O}_D$. Write $p=\pi_1\cdots\pi_k$, with $k\ge2$ and all the $\pi_i$ irreducible. Taking norms, $p^2=N(\pi_1)\cdots N(\pi_k)$, and so $k=2$ and $N(\pi_1) = N(\pi_2) = \pm p$. By Lemma \ref{lem:primenorm}, both $\pi_1$ and $\pi_2$ are prime, and so $p$ factors into primes after all, an absurdity.

Let \[ m(X) = X^2 - \sigma X + \frac{\sigma-D}{4} \in \mathbb{Z}[X] \]
be the minimal polynomial of $\tau$. Then $\mathbb{Z}[\tau]\cong \mathbb{Z}[X]/(m(X))$
 and $\mathbb{Z}[\tau]/(p) \cong \mathbb{F}_p[X]/(m(X))$. Since $p$ is not prime in $\mathcal{O}_D$, the quotient ring $\mathbb{Z}[\tau]/(p)$ is not a field, and so $m(X)$ factors nontrivially over $\mathbb{F}_p$. Thus, for some integers $x$ and $x'$,
 \begin{equation}\label{eq:mfact} m(X) \equiv (X-x) (X-x') \pmod{p}. \end{equation}
 Comparing coefficients of $X$ on both sides, we find that $x+x' \equiv \sigma\pmod{p}$, and so we can assume that
 \[ \frac{p+\sigma}{2} \le x\le p \text{\quad if $D>0$}, \qquad\text{and}\qquad \sigma\le x \le \frac{p+\sigma}{2} \text{\quad if $D<0$}. \]
By \eqref{eq:mfact}, $m(x) \equiv 0\pmod{p}$. Moreover, our inequalities for  $x$ guarantee that
\[ |m(x)| < p^2. \]
Indeed, if $D>0$, then (keeping in mind \eqref{eq:plower})
\[ p^2 > m(x) = \left(x-\frac{\sigma}{2}\right)^2 - \frac{D}{4} \ge \frac{p^2-D}{4} > -p^2, \]
while if $D < 0$, then
\[ 0 < m(x) = \left(x-\frac{\sigma}{2}\right)^2 + \frac{|D|}{4} \le \frac{p^2+|D|}{4} < p^2. \]
Write $m(x)=pr$, where $|r| < p$. By the minimality of $p$, every prime dividing $r$ factors into primes of $\mathcal{O}_D$, and so $r$ itself factors, up to sign, as a product of primes of $\mathcal{O}_D$. Thus, for some primes $\eta_1,\dots,\eta_\ell$ of $\mathcal{O}_D$,
\[ (x-\tau)(x-\bar{\tau}) = m(x) = \pm p \eta_1\cdots\eta_\ell. \]
Since $\eta_1$ is prime, $\eta_1$ divides either $x-\tau$ or $x-\bar{\tau}$. Divide both sides of the equation by $\eta_1$ and continue the process with $\eta_2$. Eventually we are led to a factorization of the form
\[ \frac{x-\tau}{\Pi_1} \cdot \frac{x-\bar{\tau}}{\Pi_2} = \pm p, \qquad\text{where}\qquad
\Pi_1 = \prod_{i \in \mathcal{I}} \eta_i,\quad \Pi_2 = \prod_{i \in \mathcal{I}^{c}} \eta_i
\]
 for some $\mathcal{I} \subset \{1,2,\dots,k\}$, where $\mathcal{I}^{c} = \{1,2,\dots,k\}\setminus \mathcal{I}$. Multiplying by $\pm 1$ if necessary, we obtain a factorization of $p$ as $\alpha\beta$, say. If $\alpha$ or $\beta$ is a unit, then the other is a unit multiple of $p$. But that implies  $p\mid x-\tau$ or $p\mid x-\bar{\tau}$, which is absurd. (Both $\{1,\tau\}$ and $\{1,\bar{\tau}\}$ are $\mathbb{Z}$-module bases of $\mathcal{O}_D$, and so when a multiple of $p$ is written as $a+b\tau$ or $a+b\bar{\tau}$, both $a$ and $b$ must be multiples of $p$.) So $\alpha,\beta$ are nonunits. Now taking norms shows that $N\alpha = N\beta=\pm p$, so that $\alpha,\beta$ are prime by Lemma \ref{lem:primenorm}. Thus, $p$ has a factorization into primes of $\mathcal{O}_D$ after all, contradicting the choice of $p$.
\end{proof}

\begin{rmk} In 1912/1913, Frobenius \cite{frobenius12} and Rabinowitsch \cite{rabinowitsch13} (independently) published the following striking result: For each integer $q \ge 2$,
\begin{multline*} x^2-x+q~\text{ is prime for all integers $0 < x < q$}\\\text{if and only if}\quad \text{$\mathbb{Z}[\tfrac{1}{2}(1+\sqrt{1-4q})]$ is a UFD}; \end{multline*}
see \cite[Chapter 11]{pollack17} for an exposition. For example, since $\mathbb{Z}[\frac{1}{2}(1+\sqrt{-163})]$ is a UFD, the polynomial $x^2-x+41$ assumes prime values for $x=1,2,\dots,40$. The  ``only if'' half of the proof is the more difficult of the two, and for this most modern treatments fall back on the theory of the class group. Theorem \ref{thm:main} allows one to fashion a completely elementary proof (apply Theorem \ref{thm:main} in place of Proposition 11.13 in \cite{pollack17}; alternatively, Ramirez V.'s Theorem 3.1 from \cite{ramirez16} can be used).  Indeed, these arguments prove a sharper version of the forward direction, which has the following consequence: $x^2-x+41$ being prime for just $x=1,2,3,4$ implies that $x^2-x+41$ must continue being prime all the way to $x=40$. Certain relatives of Rabinowitsch's theorem for real quadratic orders can  be given elementary proofs in a parallel way (compare with \cite{ramirez19}).
\end{rmk}

\subsection*{Acknowledgements}
The authors are supported by the National Science Foundation (NSF) under awards DMS-2001581 (P.\,P.) and DMS-1454767/DMS-2000093 (N.\,S.). They thank Enrique Trevi\~no and the referees for helpful suggestions. In particular, they are grateful to a referee for pointing out that the argument applies for orders other than the maximal one.