\section{Introduction}
Up-scattering of low-energy photons by Inverse Compton effect on a hot gas of
electrons (i.e. Comptonization) is a common mechanism for high-energetic
plasmas, especially for accreting sources such as Black Hole Binaries (BHBs) and Active
Galactic Nuclei (AGN). It is widely believed to be
the main source of the hard X-ray emission we observe in radio quiet sources (i.e. AGN without a jet): soft photons produced by the accretion disc are
Comptonized by a medium of hot electrons (usually referred to as the
``X--ray corona") whose geometry and physical parameters are mostly unknown.
If the electrons are thermally distributed in the corona with a certain energy,
$kT_e$, the resulting X-ray spectrum produced by Comptonization, usually
referred to as the primary component, will be a power-law with a cut-off
proportional to the thermal energy of the corona depending on the optical depth
of the corona, $\tau$. For some specific geometry and optical thickness of the
corona, this proportionality has been found to be $E_{cut} \sim 2-3$ $kT_e$ (e.g. \cite{POP2001}).

The primary component irradiates the inner accretion disc and produces
a reflection spectrum whose main features are the fluorescent Iron
K$\alpha$ at 6.4--6.7 keV (depending on the ionization of the disc) and
a 'Compton hump' due to reflection of $\gtrsim 10$ keV photons from the disc.
Since reflection may originate from disc parts which are close to the black
hole (BH), general relativistic (GR) effects (most prominently the gravitational
redshift) may affect the observed spectrum. The GR effects on the reflection
X-ray spectral components have been extensively used for spin and corona
size measurements (for redshifted iron line e.g. \cite{Fabian2002}, \cite{Miniutti2007}, and for redshifted Compton hump e.g. \cite{Svoboda2015}).
Most of these measurements have suggested the presence of highly spinning black
holes and compact coronae in AGN.  X-ray spectral-timing and reverberation
analysis of AGN spectra support this conclusion and in many cases the position
of the corona is thought to be less than $3-10 r_g$ above the black hole
(e.g. \cite{Kara2016},  \cite{Emmanoulopoulos2014}, \cite{DeMarco2011}).

In the last years, thanks to the large broadband observations obtained mainly
by \textit{NuStar}, spectra up to hundreds of keV have been collected for
many AGN and BHBs. Some of them show a clear high energy cut-off in
the power-law, allowing to estimate the thermal energy of the corona, $kT_e$ (for an exhaustive collection of results consult \cite{Fabian2015} and \cite{Lubinski2016}).  If the X--ray source lies close to the BH, GR effects such as gravitational redshift, relativistic beaming (if the corona rotates) and light bending will become important and may affect the estimation of thermal energy of the corona derived from the observed high-energy cut-off.  GR effects to the X--ray primary component have already been studied in the case of point-like sources located on the axis of disc rotation by  \cite{Niedzwiecki2016}, \cite{Fabian2015} and \cite{Fabian2017}.

Point-like and extended disc-like coronae orbiting or rotating above the disc (e.g. \cite{Wilkins2012} and  \cite{Wilkins2016}) have been extensively studied in the last few years and GR corrections are usually taken into account for such geometries in order to properly calculate emissivity profiles from the disc reflecting the X-ray primary emission coming from such coronae as well as to properly model the expected X-ray reverberation signal. However, in many works involving spectral analysis only, GR corrections on primary component are neglected.

The main goal of this work is to provide correction factors for the observed energy cut-offs in the case of X--ray spectra emitted by  a disc-like corona, i.e.  a flat, extended corona over the surface of the accretion disc. The spectrum emitted is affected mainly by gravitational redshift and Doppler shifts. We compute and provide such factors for a large number of coronal radii, inclinations, constant and radially decreasing emissivity profiles, for a maximally rotating and a Schwarzschild black hole. Depending on the inclination and the size of the corona, the intrinsic energy cut-offs can be significantly larger, but also {\it smaller}, than the observed ones. They can also be significantly larger than the respective correction factors in the case of point-like sources. Our results can be used in the study of the intrinsic energy cut-off (and hence coronal temperature and intrinsic thermal energy) distribution in AGN by simply apply the corrections we furnish to the observed energy cut-off derived by usual spectroscopical analysis.
Comparison of the resulting intrinsic distributions may offer additional clues regarding the geometry of the X-ray sources in these objects.

\section{Relativistic effects to the continuum spectrum}
We assumed two compact, simple geometries for the X--ray emitting region.
The first is the so-called lamp-post geometry where the source is static and
located on the BH spin axis, at a height $h$ above the disc plane (\cite{Matt1991}).
The second is a geometrically thin corona extended
from the radius of the innermost stable circular orbit (ISCO) to an outer radius, $r_{out}$,
in the equatorial plane just above the accretion disc and rotating with
Keplerian velocity around the BH (hereafter the
``disc corona''). Both height and radii are measured in units of the
gravitational radius, $r_g=GM/c^2$. While both these geometries are highly
simplified, one can imagine the
geometry of the real X--ray compact source to be in between these two extreme
configurations (with $h$ and $r_{out}$ being representative
of the source's height and extension, respectively).

To study the GR effects on the X--ray continuum spectrum, we assumed that the X-ray source emits isotropically (in the rest-frame) with an intrinsic spectrum of the form:

\begin{equation}
\label{intr-pl}
F(E)=NE^{-\Gamma}\exp(\dfrac{-E}{E^{i}_{cut}}).
\end{equation}

Then, we computed the spectrum that a distant observer would detect, by taking
into account all relativistic effects for both geometries (gravitational and Doppler energy shifts and light bending (\cite{Karas2006}). In both cases
the spectrum will be shifted in energy and its normalisation will change while the shape of the power-law spectrum, i.e. the power-law index, $\Gamma$, will remain the same. We are interested in particular in the energy shift of $E_{cut}$, since this observable is used to measure the thermal energy of the corona. We therefore define the $g$-factor\footnote{Note that the g-factor in GR is usually defined conversely but we find this definition more useful in our case. If g-factor is known, we can multiply the observed cut-off energy by this number to get the intrinsic cut-off energy.} as the ratio of the intrinsic, $E^{i}_{cut}$, to the observed cut-off energy, $E^{o}_{cut}$, e.g.
\begin{equation}
\label{g-fact}
g = E^{i}_{cut}/E^{o}_{cut}.
\end{equation}

In the case of the lamp-post corona $g$ is independent of inclination (as only gravitational redshift contributes ) and therefore it can be calculated using the simple relation,

\begin{equation}
\label{g_LP}
g =
    \sqrt{\dfrac{h^2+a^2}{h^2-2h+a^2}},
\end{equation}
where $h$ is the height of the corona above the BH, and $a$ is the non-dimensional  spin parameter. This is a well known formula, and has already been used in the past, for example, by \cite{Fabian2015} and \cite{Fabian2017} to correct the observed cut-off energies in their sample.

In the case of a disc corona, we expect the $g-$factor to depend on the radius of the corona, its emissivity profile, as well as the inclination. However, due to the fact that the primary power-law cut-off emission from different parts of the disc corona will be shifted by different $g-$factors due to the fact that both the gravitational redshift as well as the Doppler shift (corona rotates) vary across the corona, the overall $g-$factor, as defined by eq.~(\ref{g-fact}), cannot be estimated analytically. We give below an example of how we can estimate $g-$factor in this case.

Let us consider a disc corona which extends from the horizon of a maximally rotating black hole (spin $a = 1$) to 2.5 $r_g$, with uniform emissivity and an intrinsic spectrum (in its rest frame) defined by eq.~(\ref{intr-pl}), for $\Gamma=2$ and $E^{i}_{cut}=100$ keV (green solid line in Fig.~\ref{spectrum}).  The red solid line in the same figure shows the integrated spectrum, as it would be detected by a distant observer with an inclination\footnote{We define inclination as the angle, $\theta$, between the line of sight and the BH spin axis} of $\theta=30 ^{\circ}$. Both spectra were artificially renormalized to the value of 1 at 1 keV, to better see the energy shift and the shape of the observed spectrum. We then fitted the observed spectrum (over the $1-100$ keV band) with a model of the form,

\begin{equation}
\label{obs-pl}
F_{obs}(E) = NE^{-\Gamma}\exp\left(g\,\dfrac{-E}{E^{i}_{cut}}\right).
\end{equation}
\noindent
{The normalisation constant, $N$, power-law index, $\Gamma$, and the energy shift of cut-off, $g$, were free parameters during the fitting process. The blue dashed line in Fig.~\ref{spectrum} shows the best-fit line. The best-fit $g-$factor is 2.08. In other words, for this disc corona, the observer will detect an energy cut-off which will be {\it half} of $E^i_{cut}$. Strictly speaking, the shape of the observed spectrum (red line) is not identical to a power-law cut-off model as defined by eq.~(\ref{obs-pl}). The drop due to the cut-off energy in the observed spectrum is slightly broadened (relativistically smeared), compared to an exponential cut-off. Nevertheless, for all practical purposes, the exponential cut-off model (blue dashed line) describes the shape of the observed spectrum very well. We therefore accept the best-fit $g-$value as representative of the $g-$factor for this disc corona.

\begin{figure}
\includegraphics[width=88.mm,clip=]{figure1.eps}
\caption{\label{spectrum} The intrinsic spectrum emitted from a disc corona
extending from the horizon of a maximally rotating black hole (spin = 1) to 2.5
$r_g$ (red line), and the spectrum detected by a distant observer  at
$30 ^{\circ}$  from the spin-axis (green line). The dashed-blue line indicates
the best-fit of the function defined by eq.~(\ref{obs-pl}) to the observed
spectrum.}
\end{figure}

\section{The observed high-energy cut-off}

In order to determine the $g-$factor in the case of disc coronae, we followed the procedure described in the previous section for a large number of $r_{out}$ and inclinations.  We considered an extreme Kerr ($a=1$) and a Schwarzschild BH ($a=0$), with an inner disc radius of 1 and 6 $r_g$, respectively. We also considered the case of a disc corona with a constant emissivity and a radially reduced emissivity, as $r^{-3}$. In all cases we assumed that the intrinsic spectrum of the corona is described by eq.~(\ref{intr-pl}). We estimated the spectrum that would be observed at infinity (the ``model-observed spectrum"), taking into account all the GR effects, and we renormalized its value to 1 at 1 keV. We fitted it with the model defined by eq.~(\ref{obs-pl}), and we accept the best-fit $g-$value ($g_{fit}$) as the $g-$factor for the disc corona, with any given $r_{out}, \theta$, and emissivity profile.

During the fits, we always let $N$ and  $\Gamma$ as free parameters but, in any case,  $N_{fit}$ and $\Gamma_{fit}$ would always be almost identical to 1 and to the intrinsic slope, $\Gamma^i$. We fitted the model-observed spectra in the energy range 1--100 keV, which is broadly similar to the energy range that people would use in practise, with combined $XMM-$Newton and $NuSTAR$ data. Increasing (up to 1 MeV) or decreasing (down to 0.2 keV) the limits of the fitting energy band does not affect the resulting $g_{fit}$. However, the best-fit $g-$value changes slightly if we move the lower energy limit from 1 to 3 keV. The difference is less than 8\% in the case of very compact coronae ($r_{out}< 1.5 \, r_g$), and only for a few inclinations (in most cases the difference is less than $\sim 2-3$\%).

Since the model-observed spectrum is not entirely identical to the intrinsic spectrum, as we discussed above, we tried various $\Gamma^i$ and $E^i_{cut}$ values for all the $r_{out}$ and $\theta$ combinations we considered. In this way we can investigate how robust the $g_{fit}$ values are to different spectral shapes. We tried the values of $E^i_{cut}=30, 100, 300$ and 1 MeV, and we found that $g_{fit}$ does not vary substantially, as long as $E^i_{cut}>100$ keV. For smaller intrinsic cut-off values, we observed a difference by no more than $\sim 8$\% (in the case of the most compact disk coronae). In the case of different intrinsic spectral values, we found that $g_{fit}$ increases by $14$\% (at most), from $\Gamma^i=3$ to $\Gamma^i=1$, when $r_{out}<1.5$.

Table \ref{table:1} lists our results in the case of a disc corona, with $\Gamma^i=2$ and $E^i_{cut}=100$ keV. Rows list the $g-$factor for different inclinations ($1^{\circ}<\theta<85^{\circ}, \Delta\theta=1^{\circ}$) and columns list the $g-$factor for various $r_{out}$ (from 1 or 6 $r_g$ up to to 20 $r_g$). Both the constant and $r^{-3}$ emissivity results are listed. The lamp-post results are given in the last row for different corona heights (equal to the $r_{out}$ listed on the top row). Their estimation is simply based on eq.~(\ref{g_LP}), but we list them for completeness. Given the differences we found for various $\Gamma^i$, we provide similar Tables for $\Gamma^i=1, 1.5, 2, 2.5, 3, 3.5$ (and  $E_{cut}^{i}=100$ keV), online.

\begin{table}
\caption{$g-$factors for disc coronae with $\Gamma^i=2, E^i_{cut}=100$ keV, and constant/radial power-law ($r^{-3}$) emissivity, for different observer inclination (rows), and corona outer radius (columns). Lamp-post (LP) $g-$factors are listed in the last row, for different heights (equal to $r_{out}$, as listed in the the first row). The full table is available on line. }
\label{table:1}
\centering
\begin{tabular}{l | c | c | c | c | c}
\hline
\hline
$r_{out}$   & 1.52 & 2.50 & 5.05  & 10.2  & 20.0  \\
$\theta$ &  &  &  &  &  \\
\hline
\vspace{2mm}
 5$^\circ$&$\nicefrac{6.00}{6.21}$&$\nicefrac{2.66}{2.88}$&$\nicefrac{1.61}{1.84}$&$\nicefrac{1.27}{1.52}$&$\nicefrac{1.14}{1.40}$\\
\vspace{2mm}
30$^\circ$&$\nicefrac{4.26}{4.41}$&$\nicefrac{2.08}{2.24}$&$\nicefrac{1.40}{1.58}$&$\nicefrac{1.18}{1.38}$&$\nicefrac{1.09}{1.31}$\\
\vspace{2mm}
60$^\circ$&$\nicefrac{2.10}{2.18}$&$\nicefrac{1.24}{1.33}$&$\nicefrac{1.04}{1.14}$&$\nicefrac{1.00}{1.09}$&$\nicefrac{1.00}{1.08}$\\
\vspace{2mm}
80$^\circ$&$\nicefrac{1.15}{1.19}$&$\nicefrac{0.88}{0.93}$&$\nicefrac{0.89}{0.91}$&$\nicefrac{0.93}{0.91}$&$\nicefrac{0.96}{0.92}$\\
\hline
LP & 3.50 & 1.80 & 1.27 & 1.11 & 1.05 \\
\hline
\hline
\end{tabular}
\end{table}

\begin{figure}
\includegraphics[width=88.mm,clip=]{figure2a.eps}
\includegraphics[width=88.mm,clip=]{figure2b.eps}

\caption{\label{g}{\it Top panel:} The $g-$factor plotted as function of $r_{out}$, for various inclinations, in the case of a disc corona with $\Gamma^i=2, E^i_{cut}=100$ keV, constant emissivity for extreme Kerr and Schwarzschild BH (solid and dashed lines, respectively) and in the case of a disc corona with a $r^{-3}$ emissivity profile (dotted lines). The black solid line shows $g-$factor in the case of lamp--post geometry (in this case, the $r_{out}$ values on $x-$axis correspond to $h$).
{\it Bottom panel:} Plot of $g-$factor as a function of the inclination for
several $r_{out}$ in the case of disc coronae with constant emissivity.}
\end{figure}

Fig.~\ref{g} shows a summary of our results (listed in Table \ref{table:1}). In the lamp-post case (solid
black line in the top panel), the $g-$factor does not depend on the inclination,
and it is always larger than one, since gravitational redshift is the main GR
effect that affects the spectrum. As a result, the energy of all spectral
features should be smaller in the observer's frame. The effects to the spectrum
become noticeable when the source is very close to the BH. When $h<3 \, r_g$,
the observed cut-off energy could be 40\% smaller than the intrinsic value.
It could be less than half, when $h<2 \, r_g$.

In the case of disc corona, the $g-$factor is determined mainly by gravitational
redshift and Doppler shift. The results for disc coronae with a $r^{-3}$ emissivity profile and with uniform emissivity are plotted with the dotted and solid lines, respectively. For the same $r_{out}$, the $g-$factors are always larger in the former case. This is due to the fact that  the effective radius of the $r^{-3}$ coronae is much smaller than $r_{out}$. The difference increases with increasing $r_{out}$.

The effect of Doppler shift dominates at large
inclinations ($\theta>60^{\circ}$; blue and brown lines in top panel of Fig.\,\ref{g}). In fact, at very high inclinations, the observed high energy cut-off could be up to $\sim 20 \%$ {\it larger} than
$E_{cut}^i$ (peaking at $r_{out} \sim 3 r_g$). On the other hand, the $g-$factor
could be as high as $2-6$ in the case of X--ray sources which are seen face on
($\theta<30^{\circ}$) and are close to the BH ($r_{out}\lesssim 2.5 r_g$; bottom panel in Fig.\,\ref{g}).

In the case of Schwarzschild BHs (dashed lines in the top panel of
Fig.~\ref{g}), the contribution to large $g-$factor  from radii smaller than $6 \, r_g$ is missing. Therefore, the overall $g-$factor is smaller than  the $g-$factor in the case of disc coronae around maximally rotating BHs, with the same radius. The effect decreases with increasing outer radius,  since the relative contribution of the inner disc decreases with increasing coronal size (mainly in the case of the constant emissivity disc coronae).

\section{Discussion}

In this paper we studied in detail how GR effects affect the observed cut-off energy (and hence temperature) measurements in AGN  the case of disc-like coronae. Many studies of extended X--ray coronae in AGN have been published the last few years. For example, \cite{Wilkins2012} proposed a disc-like corona in 1H 0707-495, which is located as low as $2 r_g$ above the plane of the accretion disc and extends outwards from the rotation axis to around $30-35 r_g$, based  on the comparison between observed emissivity profiles to those computed theoretically for different locations and geometries of the X--ray source, and on the study of the time lags spectrum of the source. This X--ray corona geometry is almost identical to the geometry we study in this work: a flat corona, located above the disc, with an inner radius of $r_{ISCO}$ (for a maximally rotating and a Schwarzschild BH) and an outer radius of $r_{out}$.

We computed $g-$factors for coronae with constant emissivity and with an emissivity that decreases with radius as $r^{-3}$, for a large number of inclinations ($1^{\rm o}<\theta<85^{\rm o}$) and corona outer radius ($1.1 r_g \le r_{out} \le 20 r_g$). Our estimates are valid for any intrinsic cut-off energy (as long as $E^i_{cut}\ge 30$ keV). We used a simple phenomenological model for the continuum spectrum (namely a simple power-law with high energy cut-off; \textit {cutoffpl} in XSPEC terminology). This is not identical to more physical Comptonization models such as \textit{comptt} or \textit{compps} (e.g. \cite{Matt2015}, \cite{Marinucci2014}). We chose \textit{cutoffpl} because it was easy and fast to compute our results. Nevertheless, the differences between this model and other Comptonization models is much smaller than the difference between \textit{cutoffpl} with $\Gamma^i=3$ and $\Gamma^i=1$ that we considered in our study. Therefore, the ratio of the intrinsic to the observed roll-over energy at high energies should always be given by the $g-$factors listed in Table \ref{table:1}, irrespective of the model that is used to determine $E^o_{cut}$ (or $kT_e^{o}$).

The $g-$fit values depend, slightly, on $\Gamma^i$. We observe a moderate difference of $\sim 14$\% in the $g-$fit values in the case when $\Gamma^i=1$ and 3 (and $r_{out}$ is smaller than 1.5 $r_g$). These are rather extreme situations, and the $g-$factors differ by a much smaller amount in all other cases. Nevertheless, we provide $g-$factors for a large number of intrinsic slopes, although the corrections we list in Table\ref{table:1} should be adequate in most practical situations.

In the case of lamp-post geometry, GR corrections to $E^i_{cut}$ are easily applicable, and have already been discussed by \cite{Fabian2015} and \cite{Fabian2017}. We provide $g-$factors in this case as well, for completeness. Our results show that the $g-$factor for low inclination disc coronae ($\theta<30^{\circ}$) is larger than the lamp-post $g-$factor (for the same $h$ and $r_{out}$), even in the case of constant emissivity (solid green and red lines in the top panel of Fig.\,\ref{g}).
This is due to the fact that the contribution to the corona emission from radii below $r_{out}$ is quite substantial. While the Doppler shift is low at these face on systems, the strong gravitational redshift effects at small radii can produce stronger cut-off energy shifts. On the other hand, $g-$factors in disc coronae are smaller than the $g-$factors for point-like sources in high inclination systems, when Doppler effects dominate.

In general, we find that the GR effects to the X--ray continuum spectrum can be quite significant. If the X--ray corona in AGN is extended over the disc and $r_{out}$ is small, $E^i_{cut}$ (and hence the temperature) has been systematically underestimated. For example, in the case of disk coronae with $\theta<40^{{\rm o}}$ (i.e. Type I AGN) and $r_{out}<2.5 r_g$, $E^i_{cut}$ will be at least 2 times larger than $E^o_{cut}$. The difference could be even larger than 6 in the case of ultra compact coronae in face-on objects. The difference between $E^i_{cut}$ and $E^o_{cut}$  is less than 1.5 when $\theta > 60^{\rm o}$ (this corresponds to Type II AGN). In fact, at inclinations larger than 70 degrees, the intrinsic cut-off energy could be {\it smaller} then $E^o_{cut}$, by up to 20\%.

So far, cut-off energies are estimated by fitting observed spectra without applying GR corrections to the continuum models. As we argued above, irrespective of the continuum model used, the $g-$factors we provide can be multiplied with the observed $E^o_{cut}$ in order to estimate the intrinsic cut-off energies (or temperatures) to estimate $E^i_{cut}$. Given the differences between the $g-$factors in the case of point-like and disc-like, extended coronae, it would be useful to apply both corrections, when comparing estimated temperatures with theoretical predictions (i.e. with the pair production limit). In both cases, it is necessary to assume the coronal size (either $h$ or $r_{out}$) and $\theta$. In fact, comparison of the resulting $E^i_{cut}$ with maximum allowed temperatures in X--ray corona, can put constrains to the geometry of the source (i.e. $r_{out}$ and/or $h$, depending on the assumed geometry), in addition to any other constrains that may have been estimated by fitting X--ray reprocessing features in the energy spectra, or by timing studies.

In addition to a source by source basis, our results could also be used in other ways to investigate the X--ray source in AGN. For example, cut-off energies have been measured in many AGN in the last 10--20 years or so. The results show a distribution of values  between $\sim 40-200$ keV (even higher in some cases). Our results show that, even if the intrinsic cut-off energy is the same in all AGN, there should be a considerable scatter in the observed energies, if the objects are observed  at various inclinations and/or the outer radius of the corona is not the same in all sources (see bottom panel in Fig. 2).  This possibility could explain some, or even all of the scatter in the observed cut-off energies in AGN.

To investigate this issue, we considered the results of \cite{Lubinski2016}. They analysed the hard X-ray spectra of 28 bright Seyfert galaxies observed with {\it INTEGRAL}, together with data from {\it XMM-Newton, Suzaku} and {\it RXTE}, and they determined the  mean coronal temperature ($kT_e$) in these objects. The solid black line in the upper panel of Fig.\,\ref{resplot} shows the distribution of the observed temperatures, $kT_e^{o}$ (Table 3 of \cite{Lubinski2016}).

\begin{figure}
\includegraphics[width=88.mm,clip=]{figure3.eps}
\caption{\label{resplot} The distribution of the observed $kT_e$ and $E_{cut}$ in AGN (black histograms in the upper and lower panels, respectively) from the
\protect\cite{Lubinski2016} and \protect\cite{Fabian2015}
$NuStar$ sample.  The red histogram in both panels shows the predicted $kT_e$ and $E_{cut}$ values that we would observe in the case when the intrinsic temperature/cut-off energy (indicated by the vertical dashed lines in the two panels) and coronal size is the same in all AGN (see text for details)
}
\end{figure}

In order to investigate whether the distribution of the observed temperatures could be consistent with the hypothesis of a single intrinsic temperature, we considered a range of intrinsic temperatures, $kT_e^{i}$, from 200 to 500 keV (with a step of 5 keV), and we used the $g-$factors listed in Table \ref{table:1} (in the disc corona case, with uniform emissivity), to predict the distribution of the observed temperatures ($kT_{e,mod}^{o}$) for every $kT_e^i$ value. We let $\theta_{min}$, $\theta_{max}$ and $r_{out}$ to be free variables in the estimation of the $kT_{e,mod}^o$ distribution, and we used the K-S test to compare each model distribution with the distribution of $kT_e^{o}$. We found that the observed temperature distribution is fully consistent ($p_{null}=0.8$) with the distribution of temperatures that we would observe  if $kT_e^i=440$ keV,  the corona is very compact in all objects ($r_{out} \le1.2 r_g$), and their inclination is uniformly distributed between $13^{\circ}<\theta<85^{\circ}$ (red histogram in the upper panel of Fig.\,\ref{resplot}). The average luminosity of the objects in the \cite{Lubinski2016} sample is $\sim 0.08$, in Eddington luminosity units (from their Table 5). Therefore, if the corona is less than 1.2 $r_g$, the average X--ray compactness in these objects, $l$,  should be larger than 1500. However,  the point $(\Theta\sim440/512=0.86 \, l=1500)$ lies above all pair lines in the $(\Theta,l)$ plane (see, e.g., Fig.\,2 in \cite{Fabian2015}).

We therefore conclude that the AGN in the \cite{Lubinski2016} sample cannot have the same intrinsic cut-off energy, as in this case the X--ray corona would be too small, and its compactness unphysically  too high.

We performed the same analysis using the observed cut-off energies of the {\it NuStar} sample of AGN in \cite{Fabian2015}. This sample is smaller than the \cite{Lubinski2016} sample, but {\it NuStar} data can provide accurate measurements of $E_{cut}$. The solid black line in the lower panel of
Fig.\,\ref{resplot} shows the distribution of $E_{cut}^i$ for this sample. The red histogram in the same panel indicates the distribution of the cut-off energies that we would observed if $E_{cut}^i=280$ keV and $r_{out}= \, 5r_g$ in all AGN, and the inclination of the objects is uniformly distributed between $17^{\circ}<\theta<77^{\circ}$. The model and the observed $E^i_{cut}$ distributions are fully consistent ($p_{null}=0.89$\footnote{Due to the small number of objects, and the weakness of the K-S test to account for differences in the wings of the distribution, the very large cut-off energy of NGC~5506 at $\sim 700$ keV does not affect the goodness of fit between the two distributions.}). Due to the higher lower limit in $\theta$, and the larger coronal size, the model distribution in this case is flatter than the model distribution plotted in the upper panel of Fig.\,\ref{resplot}. Assuming a factor of 2 to convert from the cut-off energy to the coronal temperature (this is the case if the optical depth is $\lesssim 1$, \cite{POP2001}) then $\Theta=0.27$ for these AGN. If the average X--ray luminosity of the objects in this sample is $\sim 0.08$ (same as before), then for a $5\, r_g$ corona, the compactness will be $l\sim 400$.

We therefore conclude that the intrinsic cut-off energy  can be $\sim 280$ keV in all AGN in the \cite{Fabian2015} sample.
The radius of the X--ray corona should be $\sim 5 r_g $ (in agreement with the results from recent spectral/timing studies in many of them), and the resulting $(\Theta,l)$ values will locate these AGN close to the pair line for a hemispherical corona (Fig. \,2 in \cite{Fabian2015}).

The analysis above is meant to be an example of how our results could be used in the study of the properties of AGN coronae.  The results are not very restrictive at the moment, but we should be able to perform a much more detailed analysis in the near future, as the number of sources with well measured high energy cut-offs should increase, thanks to {\it NuStar}, but also {\it Swift/BAT},  observations.

\section*{Acknowledgements}
We acknowledge financial support from the European
Union Seventh Framework Programme (FP7/2007-2013) under grant agreement n.312789.
MD would like to thank for the support from the project 17-02430S funded by Grant Agency of the Czech Republic.
JS acknowledges financial support from grants LTAUSA17095 and Inter-Inform LTI17.

\bibliographystyle{mnras}


\section{\label{sec:level1}Introduction}
\indent Neural network model chemistries (NNMCs) greatly reduce the computational effort needed to simulate chemical systems with ab-initio accuracy\cite{Snyder:2013aa,brockherde2017bypassing,snyder2012finding,li2016understanding,li2016pure,vu2015understanding,handley2010potential,chmiela2017machine,behler2011neural,behler2007generalized,shakouri2017accurate,behler2017first, han2017deep,yao2017many,yao2016kinetic,yao2017intrinsic,khaliullin2011nucleation,bartok2010gaussian, mones2016exploration,gastegger2017machine,kobayashi2017neural, carpenter2017empirical,kolb2017discovering,kruglov2017energy, lubbers2017hierarchical,mills2017deep,wu2017internal,khorshidi2016amp,shao2016communication,zhang2014effects, li2015permutationally2,medders2015representation,medders2013critical,reddy2016accuracy,riera2017toward,moberg2017molecular, conte2015permutationally,manzhos2014neural, manzhos2009fitting,malshe2010input,peterson2016acceleration, khorshidi2016amp,piquemal2017preface,cubuk2017representations, MBCoarseCsanyi,FFparafit,li2017machine}. They can be used to predict molecular properties \cite{CM, BoBs,lopez2014modeling, pilania2013accelerating, schutt2014represent, ma2015machine, nelson2012nonadiabatic, janet2017predicting, janet2017resolving, yao2017intrinsic, hase2017machine, mcgibbon2017improving, bereau2017non, grisafi2017symmetry, isayev2017universal, ghiringhelli2017learning, ouyang2017sisso, faber2017prediction, schutt2017quantum, MSLSTM, timoshenko2017supervised}, and design new materials or drugs\cite{li2017learning, ramsundar2017multitask, hachmann2011harvard,hachmann2014lead, isayev2015materials, kim2017machine, segler2017towards,olivares2011accelerated,guimaraes2017objective,wei2016neural, gomez2016automatic,jinnouchi2017predicting,ulissi2017machine,sun2017machine}. In spite of growing popularity, most neural network methods are still only used by their developers and are customized for a single application. The paucity of literature describing transferable accuracy, confusion about what physics can be reproduced, and dearth of open software is slowing adoption. This paper develops an open-source, transferable neural network model chemistry. We show that NNMCs are easily hybridized with physical contributions to molecular energies familiar from Molecular Mechanics and corrections to Density Functional Theory(DFT)\cite{grimme2006semiempirical}. This approach combines the best of both worlds, yielding predictable reproduction of physical long-range forces, but also featuring a linear-scaling inductive charge model which is cheaper than a Thole model\cite{thole1981molecular} yet more accurate than fixed charges. \\
\indent     Our group is one of several who have been pursuing transferable and black-box neural network model chemistries\cite{yao2017intrinsic, brockherde2017bypassing, smith2017ani, medders2013critical, behler2007generalized, han2017deep}. The field is growing so rapidly, that most non-practitioners cannot keep up with the capabilities of existing models and the outstanding problems. Readers may not appreciate that a model can achieve chemical accuracy for energies but have uselessly noisy forces. Models which provide energies at equilibrium, and those treating a fixed molecule or stoichiometry are now reliably produced\cite{han2017deep}. We will show that TensorMol-0.1 yields usefully accurate predictions of forces out-of-equilibrium by showing reproduction of infrared spectra which closely approximate our source model chemistry (wB97X-D, 6-311G**)\cite{chai2008long}, and molecular dynamics. We outline several tricks which are required to ensure the stability of long-time molecular dynamics.\\
\indent     Another distinguishing feature of our approach is the open-source nature of our package, which uses the TensorFlow tensor algebra system to compute descriptors and forces. Rather than a monolithic black-box TensorMol-0.1 is a modular collection of differentiable chemical models written in concise Python code. The components are easily joined together and extended. The methodology can be used to propagate dynamics for large molecules ($10^5$ atoms) with GPU acceleration on simple laptop computers. No significant expertise, force field refinement, or other interventions are needed to apply the method to a molecule of interest, so long as the elements are supported. The package is also interfaced with the I-PI path integral package\cite{ceriotti2014pi}, to allow for quantum simulations and enhanced sampling. \\

\section{\label{sec:level1}Methods}

\begin{figure}[t]
\includegraphics[width=0.45\textwidth]{scheme.pdf}
\caption{The schematic graph of TensorMol-0.1. Each element has its own charge network and energy network. The charge network predicts atomic charges that yield the ab-initio dipole moment. The Behler-Parinello type energy network produces a short-range embedded atomic energy, which is summed with the electrostatic energy and Van-Der Waals energy to predict the total atomization energy of molecules at and away from equilibrium. The whole computation is included in the TensorFlow graph. Analytical GPU-accelerated forces are symbolically generated by a single line of code.}
\label{Fig1}
\end{figure}

\indent     The community of Neural-Network model chemistry developers is rapidly improving the accuracy and generality of these reactive force fields\cite{behler2017first,  shakouri2017accurate, bartok2010gaussian, deringer2017machine, li2015permutationally2, smith2017ani}. The model of this paper includes several components which were the subject of recent developments in other groups\cite{behler2007generalized,behler2017first,morawietz2013density, smith2017ani, gastegger2017machine}. We will describe the details here from the bottom up citing prior art. Our notational convention will be that $i,j,k...$ are indices of atoms, $q_i$ is the charge on atom i, $z,x,y$ are atomic numbers, $\mathcal{A},\mathcal{B}, \mathcal{C}$ are molecules, and $\alpha, \beta... $ are indices of basis functions which are a product of radial and angular functions. If a function depends on all the atomic coordinates of a molecule it will be written as a vector, and those which depend on only few will be given explicit indices. The energy of TensorMol-0.1 is expressed as a sum of a short-range embedded N-body potential\cite{behler2017first}, and long-range electrostatic potential and Van-Der-Waals force:
\begin{eqnarray}
E(\vec{R}) =&&\sum_i E^\text{BP}_{z_i} (S_\alpha(\vec{R})) +\sum_{ij} E^\text{DSF}_{ij}(q_{ij}(S_\alpha(\vec{R})),R_i,R_j)\nonumber\\&&+E^\text{VDW}(\vec{R}_{ij})
\end{eqnarray}
In the above expression $E_{z_i}$ is a Behler-Parinello type energy network for the element $z$ for atom $i$. This n-body potential takes as its argument $S_\alpha$ the modified symmetry functions of Isayev and coworkers\cite{smith2017ani}:
\begin{eqnarray}
S_\alpha (radial) = && \sum_{j\neq i}e^{{-\eta(R_{ij}-R_{s})^2}}f_c(R_{ij}) \\
S_\alpha (angular) =&& 2^{1-\zeta }\sum_{j\neq i, j\neq k}(1+cos(\theta_{ijk}-\theta_s))^{\zeta}\nonumber\\&&\times e^{-\eta(\frac{R_{ij}+R_{ik}}{2}-R_s)^2}f_c(R_{ij})f_c(R_{ik})
\end{eqnarray}
Modern machine learning frameworks provide automatic differentiation of tensor algebraic expressions, allowing a force-field developer to obtain the gradient of a molecular potential $\frac{dE(\vec{R})}{d\vec{R}}$ in a single line of code, once the expression for $E(\vec{R})$ has been written. An important feature of our code is that this symmetry function is coded within the TensorFlow system\cite{tensorflow2015-whitepaper}, so all the parameters of this descriptor can be variationally optimized alongside the network weights. Our implementation of the symmetry function employs a list of nearest-pairs and triples within radial cutoffs such that the scaling of the overall network is asymptotically linear. On an ordinary laptop equipped with only a CPU a force/energy call on 20,000 atoms is less than a minute. \\
\indent     The second term of our energy expression is the Damped-Shifted Force (DSF) Coulomb energy of Gezelter and coworkers\cite{fennell2006ewald}. The charges are obtained from a sub-network which reproduces molecular dipole moments. Our charge model enforces conservation of total charge by evenly spreading any required neutralizing charge over the entire molecule or unit cell. The Damped-shifted force ensures long range continuity and differentiability of the effective Coulomb potential with smooth cutoffs. We modify the DSF kernel at short range with an "elu" type non-linearity, such that the forces within the radius of the Behler-Parinello symmetry function smoothly approach zero avoiding singularities and interference with the Behler-Parinello many-body potential.
\begin{eqnarray}
E_{DSF}= &&\left\{\begin{matrix}  E_\text{DSF(Original)}  \hfill   R > R_\text{switch}
\\ q_iq_j(a_\text{elu}e^{R-R_\text{switch}}+\beta_\text{elu})\hfill R< R_\text{switch}
\end{matrix}\right.
\end{eqnarray}
where $E_\text{DSF(Original)}$ is the energy of DSF kernel\cite{fennell2006ewald},  $R_{switch}$ is the short range cutoff for the "elu" kernel. $\alpha_\text{elu}$ and $\beta_\text{elu}$ are chosen so that the value and the gradient of $E_\text{DSF}$ are continuous at $R_\text{switch}$.
The damped-shifted force is well-suited to combination with Neural Network models because it requires no Fourier transformation to treat periodic systems with linear scaling and maps well onto TensorFlow. The last term is the van-der waals energy, which is calculated by following Grimme's C6 scheme \cite{grimme2006semiempirical}.\\
\indent We employed a two step training approach. First, the charge networks are trained to learn the atom charges that predict the dipole moment. The loss function can be written as following:
\begin{align}
 L_{dipole} = &&\sum_\mathcal{A} (\frac{\mu^\text{DFT}_\mathcal{A} - \mu^{NN}_\mathcal{A}(q_i, q_j,...)}{N_\text{atom}})^2
\end{align}
After the charge training is converged, we train the energy network. During the energy network training, the weights in charge networks are kept frozen, but they are still evaluated to calculate the electrostatic energy that is added to construct the total energy. Our Behler-Parinello many-body potential also absorbs the shape of the transition between many-body and electrostatic regions. The learning target for the energy network includes both the DFT energy and DFT force. The loss function for the energy network training is:
\begin{align}
 L_\text{energy} =&& \sum_\mathcal{A} (\frac{E^\text{DFT}_\mathcal{A} - E^\text{NN}_\mathcal{A}}{N_{atom}})^2 + \gamma \sum_\mathcal{A} (\frac{F^\text{DFT}_\mathcal{A} - F^\text{NN}_\mathcal{A}}{N_\text{atom}})^2
\end{align}
where $E^\text{NN}$ is obtained according to equation 1, $F^\text{NN}$ is calculated by taking the gradient of $E^\text{NN}$ respect to the coordinates of the atoms. $N_\text{atom}$ is the number of the atoms in the system and $\gamma$ is a parameter that controlling the portion of force loss. We employ $\gamma=0.05$.
\indent  We trained two neural networks based on two sets of data. One network ("water network") was trained on a database that includes $\sim$370,000 water clusters with 1 water molecule to 21 water molecules. The other network was trained on $\sim$3,000,000 different geometries of $\sim$ 15,000 different molecules that only contains C, H, O and N and up to 35 atoms. Since these 15K molecules were sampled randomly from the chemspider database, therefore we will refer this network as "chemspider network" in the following text. The training data were sampled using metadynamics and calculated by Qchem package\cite{shao2015advances} with WB97X-D\cite{chai2008long},  exchange correlation functional and 6-311G** basis set.\\

Each charge network and energy network contains three fully-connected hidden layers with 500 hidden neurons in each layer. For chemspider network, a network with three hidden layers with 2000 hidden neurons in each layers is used for each charge network and energy network. L2 regularization and dropout \cite{srivastava2014dropout} on last layer were used to prevent overfitting with a dropout probability of 0.3. We chose a softplus function as the non-linear activation function after extensive experimentation, and used the Adaptive moment quasi-Newton solver (Adam)\cite{kingma2014adam} to fix the weights of the network. The test sets were separated from training data by choosing a random 20\% of molecules at the outset which were kept independent throughout. Besides water we will present calculations from molecules strictly absent from either the training or test set.

\indent     To obtain scalable efficiency TensorMol uses neighborlists within cutoffs. This allows double precision energy, charge, force calculations of up to 24,000 atoms to execute in less than 90 seconds on a 2015 Intel i7 2.5GHz MacBook pro (Fig. \ref{fig:timings}). Periodic evaluations are achieved by tessellation of a unit cell with summation of energies for atoms within the cell. This results in  roughly a factor of three overhead in computational time. Speedups are obtained automatically for computers with GPUs or single-precision calculations.
\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{timing_water.pdf}
    \caption{Aperiodic timings of an energy, charge, force call for cubic water clusters at a density of 1 gm/cm$^3$. The largest $\sim60$ Angstrom cube is 4x larger than the electrostatic cutoff. The slope of a log-log version of this curve is near unity, indicating the wall-time scaling of TensorMol.}
    \label{fig:timings}
\end{figure}

\section{\label{sec:level1}Results}
The root mean square error (RMSE) on the independent test set of the energy is 0.054 kcal/mol per atom and the RMSE of the force is 0.49 kcal/mol/\AA. The left panel of figure \ref{fig:trimer} plots the potential energy surface (PES) of a water trimer when one of the water is pulled away from the other two. One can see our neural network PES is not only in good agreement with the PES of target method but also smooth. To achieve this we use a variation of the soft-plus neuron rather than the rectified linear units which are popular in computer science. The latter train more efficiently, but produce discontinuous forces. \\
\indent     The right panel shows the fraction of each of the three energy components in equation 1 to the binding energy along the trimer dissociation coordinate.  At short range, most of the binding energy is contributed by the N-body neural network potential. When the distance between the monomer and the dimer approach the cutoff distance of the neural network, the contribution of neural network potential starts to decrease and the contribution of electrostatic potential increases. After 6\ \AA\ where the neural network symmetry functions on the atoms in the monomer have no contribution from the dimer, the neural network force drops smoothly to zero and the electrostatic interaction dominates. The small difference in the energy at 7\ \AA\ is due to the difference between the Madelung energy given by the learned charges, and the genuine physical cohesive force at this distance. The dimer and monomer are beyond the symmetry function sensory radius, and so the charges are constant in this region. Future iterations of the charge network will use local-field information to improve this region of the PES. The learned inductive charges are of high quality considering their linear scaling cost. Figure \ref{fig:dimer} shows the PES and dipole change of a water dimer when the hydrogen bond is broken by rotating the OH bond. Both the PES and dipole change fit well with the DFT results. \\

\begin{figure}
    \centering
    \includegraphics[width=0.42\textwidth]{trimerbinding.pdf}
    \caption{Left panel: PES of water trimer when one water is pulled away from the other two. Right panel: Percentage contribution of binding energy between the water that is pulled away and the other two water from Behler-Parrinello atom-wise energy, electrostatic energy and van-der waals energy. Behler-Parrinello atom-wise energy contribute to most of the binding energy at the short range and electrostatic energy is the dominant contribution at long range.}
    \label{fig:trimer}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.42\textwidth]{dimerdipole.pdf}
    \caption{Top left panel: PES of breaking a hydrogen bond between two water by rotating one water around O-H bond. Top right, bottom left and bottom right panels: change of x, y, z, dipole component during the rotation, respectively}
    \label{fig:dimer}
\end{figure}

Given the increased dimension of the Hessian, it is naturally a more stringent test to reproduce forces and infrared spectra than it is to simply produce energies. The left panel and right panel of figure \ref{fig:waterir} show the optimized geometries and IR spectra of a 10 water cluster and 20 water cluster generated with our force field and DFT, respectively. Each method uses its own equilibrium geometry, so this also tests TensorMol-0.1's reproduction of non-covalent geometry. The RMSE of the distance matrix between DFT optimized geometry and TensorMol optimized geometry are 0.062\ \AA\ for the 10 water cluster and 0.180\ \AA\ for the 20 water cluster. Our force field quantitatively reproduces the DFT IR both in terms of frequencies and intensities, especially for the water bend modes and inter-monomer modes. The Mean Absolute Error (MAE) of frequencies at that those two regions are 33.2 cm\textsuperscript{-1} for the 10 water cluster and 16.2 cm\textsuperscript{-1} for the 20 water cluster. The error is slightly larger at  water OH stretching region with a MAE of 34.2 cm\textsuperscript{-1} and 13.1 cm\textsuperscript{-1}, respectively. This accuracy is comparable to high quality polarizable water force fields\cite{medders2015representation}.\\

\begin{table}
\caption{\label{tab:table1} Training details and test RMSE of each learning target. The unit of energy RMSE, gradient RMSE and dipole RMSE is kcal/mol per atom, kcal/mol/\AA\ per atom and Debye per atom, respectively.}
\begin{ruledtabular}
\begin{tabular}{lcc}
&Water Network&Chemsipider Network\\
\hline
Num of training case & 370844 & 2979162\\
Training time (days)\footnote{Training was done on single Nvidia K40 GPU}& 3  & 10 \\
Energy RMSE& 0.054 & 0.24\\
Gradient RMSE& 0.49 & 2.4\\
Dipole RMSE& 0.0082  & 0.024\\
\end{tabular}
\end{ruledtabular}
\end{table}

\begin{figure}
    \centering
    \includegraphics[width=0.42\textwidth]{water_IR.pdf}
    \caption{Simulated harmonic IR spectrum of 10 water cluster (left panel) and 20 water cluster (right panel) generated by WB97X-D/6-311G**(dashed orange line) and TensorMol force field (solid blue line).}
    \label{fig:waterir}
\end{figure}
Compared with traditional force fields, one major advantage of TensorMol is its reactivity. TensorMol is able to simulate a concerted proton transfer in a water hexamer, finding a minimum energy transition path. The PES's calculated by nudged elastic band (NEB) method \cite{henkelman2000climbing} with the TensorMol force field and DFT are shown in figure \ref{fig:waterneb}. The barrier height predicted by TensorMol is 36.3 kcal/mol, which is 6.7 kcal/mol lower than the prediction of DFT, which is remarkable considering the dearth of transition structures in the training data. Our sampling of molecular geometries uses a meta-dynamics procedure described elsewhere, so these proton transfers occur in the training data although extremely infrequently.

Encouraged by our water results, we developed a force field with applicability across the chemical space spanned by CNOH. The Chemspider dataset that we used to train our force field covers a vast chemical space containing 15 thousand different molecules and 3 millions geometries. The geometries are generated using a meta-dynamics procedure\cite{barducci2008well}, which ensures that each new geometry is a fresh part of chemical space, energies up to 400$k_bT$ are sampled in the data. We describe the details of this meta-dynamics sampling algorithm, which we have found vital to achieving robust and transferrable force-fields elsewhere. The diversity of structures makes learning the chemspider dataset a much harder task for neural networks, the test set RMSE of energy is 0.24 kcal/mol per atom and RMSE of force is 2.4 kcal/mol per atom. More importantly, the model usefully reproduces several elements of molecular structure at and away from equilibrium for molecules outside its training set. It robustly optimizes the geometries of typical organic molecules to structures that match DFT well, and yields infrared frequencies and intensities in good agreement with ab-initio calculations. It is a black-box method which does not rely on any atom type, connectivity, etc as one would need to specify in a traditional classical force-field. The few proteins we have examined remain stable and near their experimental structures when optimized or propagated at room temperature using the TensorMol-0.1 force field.\\
\begin{figure}
    \centering
    \includegraphics[width=0.42\textwidth]{water_neb.pdf}
    \caption{Reaction energy profile converged from a Nudged elastic band along the reaction coordinate of conservative proton transfer in a water hexamer cluster.}
    \label{fig:waterneb}
\end{figure}
\indent Morphine is not included in our training set. The left panel of figure \ref{fig:morphine} shows the geometry of morphine that is optimized with our force field. The RMSE of bond lengths predicted by our forcefield is 0.0067\ \AA\ and the RMSE of angles is 1.04 degrees compared with the source DFT model chemistry. The right panel plots the harmonic IR spectra generated by each method at their respective optimized geometries. One can see IR spectrum from our force field is in good agreement with the DFT IR spectrum. The MAE of our force field frequencies is 13.7 cm\textsuperscript{-1} compared with DFT frequencies. Figure \ref{fig:chemspider} shows comparisons of IR spectrum that are generated from these two methods for aspirin, typrosine, caffeine and cholesterol. All these four molecules are not included in the training set. The MAE of the frequencies predicted by our field is less than 20 cm\textsuperscript{-1} for all the four molecules compared with target DFT frequencies. The concept of a chemical bond and force constant are not enforced in any way, yet good agreement with DFT is obtained at a tiny fraction of the original cost. \\
\indent    Traditional harmonic vibrational spectra require quadratic computational effort, which works against the speed advantage of a NNMC. For large systems one can use the molecular dynamics functionality of TensorMol to simulate infrared spectra, Fourier transforming the dipole-dipole correlation function of conservative Newtonian dynamics whose cost grows linearly with the size of the system. The lower left panel of Figure \ref{fig:morphine} shows the same infrared spectrum produced by propagation in TensorMol-0.1, also showcasing the good energy conservation of TensorMol. Unlike a traditional force-field it's non-trivial to obtain smoothly differentiable NNMC's. 64-bit precision needs to be used the network cannot be made too flexible and smooth versions of typical rectified linear units need to be used. Our package can be used in this way to simulate IR of large systems with linear cost.\\
\indent     TensorMol-0.1 has a relatively simple treatment of electrostatic and Van-Der-Waals forces which we would like to augment in the future with a many-body dispersion scheme\cite{tkatchenko2012accurate}. However a main adantage of TensorMol-0.1's approach is its very low cost. No self-consistent polarization equation is solved even though the charges are inductive, and so it is easy to inexpensively calculate the electrostatic energies of even very large molecules. At shorter ranges, non-covalent interactions like hydrogen bonds are dealt with by the Behler-Parinello portion of the network. The Chemspider training data includes some examples of dimers and intra-molecular hydrogen bonds. To our surprise the treatment of inter-molecular interactions which were not targets for TensorMol-0.1 are satisfactory. Figure \ref{fig:DNA} shows the optimized geometries and binding energies of two DNA base pairs calculated by our force field. The target DFT method predicts a binding energy of 18.3 kcal/mol for the thymine-adenine (TA) pair and a binding energy of 32.4 kcal/mol for the guanine-cytosine (GC) pair. The prediction of our force field is 1.2 kcal/mol less for TA pair and 2.0 kcal/mol larger for GC pairs relative to DFT. \\
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{morphine_IR.pdf}
    \caption{Morphine geometry that is optimized by TensorMol-0.1 (upper right panel) and its harmonic IR spectrum simulated by WB97X-D/6-311G**(dashed orange line) and TensorMol force field (solid blue line) (upper left panel). Lower panels show TensorMol's real-time IR spectrum vs. DFT (left) and the conservation of energy maintained by the smoothness of the energy (right).}
    \label{fig:morphine}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{chemspider_IR.pdf}
    \caption{Harmonic IR spectrum of four different molecules simulated by WB97X-D/6-311G**(dashed orange line) and TensorMol-0.1. All the molecules are not included in the training set.}
    \label{fig:chemspider}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{DNA.pdf}
    \caption{Binding energy between the DNA base pairs vs. $\omega$B97x-D with methods at their optimized geometries. The difference between DFT and TensorMol binding energy is < 2 kcal/mol.}
    \label{fig:DNA}
\end{figure}
\indent     One holy grail of the field of Neural Network model chemistries is to simulate biological chemistry without QM-MM or bespoke force-fields. Protein simulation also demonstrates several important features of a neural network model chemistry: reasonable inter-molecular forces, stability, scalability and generalization far from small-molecule training data. TensorMol-0.1 was not trained on any peptide polymers and includes no biological data of any sort. To our pleasant surprise, even this first iteration of Neural Network model chemistry is accurate enough to perform rudimentary studies of small proteins. Figure \ref{fig:protein} shows geometries sampled from a 1 picosecond, periodic, 300K TensorMol dynamics NVT trajectory in explicit solvent. The initial structure (included in the supplement) was generated from the PDB structure 2MZX using OpenMM's automatic solvation and hydrogenation scripts\cite{OpenMM}, but includes nothing but atom coordinates. This short alpha-helix is stable, both in optimizations and dynamics, and the structures sampled during the dynamics superficially resemble the solution NMR structure. Traditional force fields will always be less expensive (by some prefactor) than NNMCs, yet the reactivity advantages of NNMCs and the ease of set up will probably lead to a rapid adoption of these methods in the biological community.
\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{protein.pdf}
    \caption{Left panel shows samples from a 1 picosecond NVT (Nos\'e) trajectory of solvated 2MZX at 300K simulated by our TensorMol force field in explicit water. Right panel is the NMR structure of 2MZX from the PDB database.}
    \label{fig:protein}
\end{figure}
\subsection{Discussion and Conclusions}
\indent     We have presented a transferable neural network model chemistry with long-range Coulombic physics, and a short-range n-Body potential, TensorMol-0.1. The model is integrated in a concise open-source python package which provides many of the types of simulation commonly used in chemistry. The method can be used to scan conformational and chemical space along the singlet neutral potential energy surface with high throughput and accuracy using nothing but atomic coordinates. TensorMol-0.1 is not the final iteration of Neural-Network model chemistry, although it shows just how easily DFT-quality predictions can be made by models with drastically lower cost. Inexpensive post-DFT corrections such as Many-Body Dispersion \cite{tkatchenko2012accurate} will become even more powerful when integrated with these potentials, opening the door to quantitative treatments of large systems. These methods will compete aggressively with DFT packages, and provide an interesting complement to QM-MM type simulations in the near future. \\
\indent     Methods like TensorMol-0.1 suggest the following directions for the future improvement of Neural-Network model chemistries:
\begin{itemize}
\item Generalize descriptors to encode other physical atom properties besides charge (spin, polarizability)
\item Develop accurate descriptors whose cost grows linearly with the number of elements treated
\item Extend the range of the N-Body embedding
\item Explore the hierarchy of physical detail between Force-Fields and semi-empirical electronic structure.
\end{itemize}
These goals must be pursed alongside honest test suites with open programs and data. NNMCs that cannot produce forces and MD trajectories with a demonstrable and compelling cost advantage over DFT should only be pursued if they offer interesting qualitative insights.

\begin{acknowledgments}
The authors gratefully acknowledge Notre Dame's College of Science for startup funding, Oak Ridge national laboratory for a grant of supercomputer resources and NVidia corporation.
\end{acknowledgments}

\input{tensormol.bbl}

\clearpage

\end{document}


\section{Introduction}
In this paper we consider a problem related to  inhomogeneous Diophantine approximation. Given $\pmb{\theta} =(\theta_1,\theta_2)\in \mathbb{R}^2 $ we study    the set of couples $(\eta_1,\eta_2)\in\mathbb{R}^2$ such that the system of two linear forms
$$
|| x\theta_1 - \eta_1|| , ||x\theta_2 - \eta_2||,
$$
where $||\cdot ||$ stands for the distance to the nearest integer
is badly approximable. We prove a statement complementary to our recent result  from
\cite{BMS}.
We construct $\pmb\theta$ such that the set
$$\mathrm{Bad}_{\pmb{\theta}} :=  \{(\eta_1,\eta_2):  \inf_{x\in\mathbb{N}} x^{\frac{1}{2}} \max_{i=1,2} \|x \theta_i-\eta_i\|>0\}
$$
is not isotropically winning.

Our paper is organized as follows. In Section 2 we discuss different games appearing in Diophantine problems. In Section 3
we give a brief survey on inhomogeneous badly approximable systems of linear forms and formulate our main result Theorem 3.1.
Sections 4 and 5 are devoted to some auxiliary observations. In Sections 6, 7, 8 we give a proof for Theorem 3.1.

\section{Schmidt's game and its generalizations}
The following game was introduced by Schmidt in \cite{Sch1}, \cite{Sch2}, \cite{Sch3}. Let $0<\alpha, \beta <1$.
Suppose that two players A and B choose by the turn a nested sequence of closed balls
$$
B_1 \supset A_1 \supset B_2 \supset A_2 \supset \ldots
$$
with the property that the diameters $|A_i| , |B_i|$ of the balls $A_i, B_i$ satisfy
$$
|A_i| = \alpha |B_i|, \qquad |B_{i+1}| = \beta |A_i| \qquad \forall i = 1,2,3,\ldots
$$
for fixed $0<\alpha, \beta <1$. A set $E\subset \mathbb{R}^n$ is called an $(\alpha,\beta)-winning$ if player A has a strategy which guarantees that intersection $\cap A_i$ meets $E$ regardless of the way B chooses to play.  A set $E\supset \mathbb{R}^n$ is called an $\alpha-winning$ set if it is $(\alpha,\beta)$-winning for all $0<\beta <1$.

There are different modifications of Schmidt's game: strong game and  absolute game introduced in \cite{Mcm}, hyperplane absolute game introduced in \cite{WK}, potential game considered  in \cite{FSU} and some others.
In our recent paper
 \cite{BMS} we introduced isotropically winning sets. Let us  describe  here some of these generalizations in more details.

The definition of an absolutely winning set was given by  McMullen  \cite{Mcm}. Consider the following game.
Suppose $A$ and $B$ choose by the turn a sequence of balls $A_i$ and $B_i$ such that the sets
$$
B_{1} \supset (B_1 \backslash A_1) \supset B_2 \supset (B_2 \backslash A_2) \supset B_3 \supset \ldots
$$
are nested.
For fixed $0<\beta <1/3$ we suppose
$$
|B_{i+1}| \ge \beta |B_i| , \qquad |A_i|\le \beta |B_i|.
$$
We say $E$ is an $absolute$ $winning$ set if for all $\beta \in (0, 1/3),$ player A
has a strategy which guarantees that $\cap B_i$ meets $E$ regardless of how B chooses to play.
Mcmullen proved that an absolute winning set is $\alpha$-winning for all $\alpha <1/2$.
Several examples of absolute winning sets were exhibited by McMullen in \cite{Mcm}.
In particular  set of badly approximable numbers in $\mathbb{R}$ is absolutely winning. However the set of simultaneously badly approximable vectors in  $\mathbb{R}^n$ for $ n > 1$
is not absolutely winning.

In \cite{BMS} another strong variant of winning property was given.
We say that a set $E \subset \mathbb{R}^n$ is {\it isotropically winning} if for each $d \le n$ and for each $d$-dimensional affine subspace $\mathcal{A} \subset \mathbb{R}^n$ the intersection $E \cap \mathcal{A}$ is $1/2$-winning
for Schmidt's game considered as a game in $\mathcal{A} $.
It is clear that an absolute winning set is isotropically winning for each $\alpha \le 1/2$.

\section{Inhomogeneous approximations}
The first important result on inhomogeneous appxohimations in  one-dimensional case is due to
Khintchine.
In \cite{KH1}
he proved that there exists an absolute constant $\gamma$ such that for every $\theta\in \mathbb{R}$ there exists $\eta\in \mathbb{R}$
such that
$$
\inf _{q\in \mathbb{Z} } q ||q\theta - \eta|| > \gamma
.$$
Later  (see \cite{KH2}, \cite{KH3}) he
 proved that for given positive numbers $n, m \in \mathbb{Z}$ there exists a positive constant $\gamma_{nm}$ such that for any $m \times n$ real matrix $\pmb{\theta}$ there exists a vector $\pmb{\eta} \in \mathbb{R}^n$ such that
$$
\inf_{\pmb{x}\in\mathbb{Z}^m \backslash \{0\}} (|| \pmb{\theta} \pmb{x} - \pmb{\eta}||_{\mathbb{Z}^n})^n ||\pmb{x}||^m >\gamma_{nm}
$$
(here $||\cdot||_{\mathbb{Z}^n}$
stands for the distance to the nearest integral point in $\sup$-norm).
 These results are presented in a wonderful book \cite{Cas} by  Cassels.

Jarnik in \cite{J1}, \cite{J2} proved a generalization of this statement.
Suppose $\psi (t)$ is a function decreasing to zero as $t \rightarrow + \infty$.
 Let $\rho (t)$ be the function inverse to the function $t \mapsto 1/\psi(t)$. Suppose that for all $t > 1$ one has $\psi_{\theta}(t) \le \psi(t)$.
Then there exists a vector $\pmb{\eta} \in \mathbb{R}^n$ such that
$$
\inf_{\pmb{x}\in\mathbb{Z}^m \backslash \{0\}} (|| \pmb{\theta} \pmb{x} - \pmb{\eta}||_{\mathbb{Z}^n})\cdot \rho(8m\cdot ||\pmb{x}||) >\gamma
$$
with appropriate $\gamma = \gamma(n,m)$.

Denote by
$$
\mathrm{Bad}_{\theta} = \left\{\alpha \in [0,1):  \inf_{q\in\mathbb{N}}q\cdot ||q \theta -\alpha||>0\right\}.
$$
It happened that
the winning property of this inhomogeneous Diophantine set was considered quite recently.
In \cite{T} Tseng showed that $\mathrm{Bad}_{\theta}$ is winning for all real numbers $\theta$ in classical Schmidt's sense.
For the corresponding multidimensional sets
$$
\mathrm{Bad}(n,m)=\left\{\pmb{\theta}\in\mathrm{Mat}_{n\times m}(\mathbb{R}):\, \inf_{q\in\mathbb{Z}^m_{\neq 0}}\max_{1\leq i\leq n}(|q|^{\frac{m}{n}}\|\pmb{\theta}_i(q)\|)>0\right\}.
$$
the winning
property  shown for example in
  \cite{ET} and \cite{NGM}.

  Further generalizations deal with the {\it twisted} sets
$$
\mathrm{Bad}(i,j)=\left\{(\theta_1,\theta_2)\in\mathbb{R}^2:\, \inf_{q\in\mathbb{N}}\max(q^i\|q\theta_1\|, q^j\|q\theta_2\|)>0\right\},
$$
where $i,j$ are real positive numbers satisfying $i+j=1$,
introduced by Schmidt.
 An in \cite{An} proved that $\mathrm{Bad}(i,j)$ is winning for the standard Schmidt game.
 In higher dimension, we fix an $n$-tuple $\textbf{k}=(k_1,\ldots,k_n)$ of real numbers satisfying
\begin{equation}\label{weight}
k_1,\ldots,k_n>0\qquad\mbox{and}\qquad \sum_{i=1}^n k_i=1,
\end{equation}
and define
$$
\mathrm{Bad}(\textbf{k},n,m)=\left\{\pmb{\theta}\in\mathrm{Mat}_{n\times m}(\mathbb{R}):\, \inf_{q\in\mathbb{Z}^m_{\neq 0}}\max_{1\leq i\leq n}(|q|^{mk_i}\|\pmb{\theta}_i(q)\|)>0\right\}.
$$
Here, $|\cdot|$ denotes the supremum norm, $\pmb{\theta}=(\pmb{\theta}_{ij})$ and $\pmb{\theta}_i(q)$ is the product of the $i$-th line of $\pmb{\theta}$ with the vector $q$, i.e.
$$
\pmb{\theta}_i(q)=\sum^m_{j=1} q_j \pmb{\theta}_{ij}.
$$
In twisted setting much less is known.
In particular up to now winning property of the set $
\mathrm{Bad}(\textbf{k},n,m)$ in dimension greater that two is not proved.

Given $\pmb{\theta}\in\mathrm{Mat}_{n\times m}(\mathbb{R})$, we define
$$
\mathrm{Bad}_{\pmb{\theta}}(\textbf{k},n,m)=\left\{x\in\mathbb{R}^n:\, \inf_{q\in\mathbb{Z}^m_{\neq 0}}\max_{1\leq i\leq n}(|q|^{mk_i}\|\pmb{\theta}_i(q)-x_i\|)>0\right\}.
$$
Harrap and Moshchevitin in \cite{HM} showed that this set is winning provided that $\pmb{\theta}\in\mathrm{Bad}(\textbf{k},n,m)$.
In \cite{BMS}
it was proved that
if we suppose that
 $\pmb{\theta}\in\mathrm{Bad}(\textbf{k},n,m)$
 the set $
\mathrm{Bad}_{\pmb{\theta}}(\textbf{k},n,m)$ is isotropically winning\footnote{In fact the approach from \cite{BMS}
gives a little bit more.
Instead of property that for any subspace $\mathcal{A}$  the intersection $E \cap \mathcal{A}$  is 1/2 winning  in $\mathcal{A}$ one can see that it is $\alpha$-winning for all $\alpha \in (0,1/2]$. It is not completely clear for the author if these two properties are equivalent. (For a closely related problem see \cite{Dremov}.)}.

We should note that even in the case
  $n=2, m=1$
  it is not known if the set
  $\mathrm{Bad}_{\pmb{\theta}}(\textbf{k},2,1)$ is $\alpha$-winning for some positive $\alpha$ without the condition
   $\pmb{\theta}\in\mathrm{Bad}(\textbf{k},2,1)$.
  .

In this article we show that the condition for $\pmb{\theta}$ be from $\mathrm{Bad}(\textbf{k},n,m)$ is essential for the isotropically winning property, and prove the following theorem.
\begin{theorem}\label{th}

There exists a vector $\pmb{\theta}= (\theta_1,\theta_2)$ such that

1) $ 1,\theta_1,\theta_2$ are lineary independent over $\mathbb{Z}$;

2) $\mathrm{Bad}_{\pmb{\theta}} :=  \{(\eta_1,\eta_2):  \inf_{x\in\mathbb{N}} x^{\frac{1}{2}} \max_{i=1,2} \|x \theta_i-\eta_i\|>0\}$ is not isotropically winning.
\end{theorem}

\section{Some more remarks}
In the sequel
$
\pmb{x} = (x_0, x_1,x_2)
$ is a vector in $\mathbb{R}^3$,
$| \cdot |$  stands  for the Euclidean norm of the vector, and by $(\bm{w}, \bm{t})$ we denote the
innerproduct of vectors $\bm{w} $ and $\bm{t}.$

\begin{comment}
\begin{theorem}\label{th}

There exists a vector $\pmb{\theta}= (\theta_1,\theta_2)$ such that

1) $ 1,\theta_1,\theta_2$ are lineary independent over $\mathbb{Z}$;

2) $\mathrm{Bad}_{\pmb{\theta}} :=  \{(\eta_1,\eta_2):  \inf_{x\in\mathbb{N}} x^{\frac{1}{2}} \max_{i=1,2} \|x \theta_i-\eta_i\|>0\}$ is not isotropically winning.
\end{theorem}
\end{comment}
The proof of the Theorem \ref{th} we will give in Section \ref{indep}.
There we will construct a special $\pmb{\theta}$ and one-dimensional affine subspace $\mathcal{P} $ such that  $\pmb{\theta} \in \mathcal{P}$ and for the segment $\mathcal{D} = \mathcal{P} \cap \{|\pmb{z} -\pmb{\theta}|\le 1\}$  one has $\mathcal{D} \cap \mathrm{Bad}_{\pmb{\theta}} = \emptyset.$
Moreover given arbitrary positive function $\omega(t)$ monotonically (slowly) increasing to infinity we can ensure that for all $\pmb{\eta} = (\eta_1,\eta_2) \in \mathcal{D}$ there exist infinitely many $x \in \mathbb{Z}$ such that
$$
\max_{i=1,2} \|x \theta_i-\eta_i\|< \dfrac{\omega(x)}{x}.
$$
To explain the construction of the proof it is useful to consider the case when $\theta_1, \theta_2,1$ are linear dependent. This case we will discuss in the Section \ref{dep}.

{\bf Remark 1.}
 From the result of the paper \cite{BMS}  it follows that the vector $\pmb{\theta}$ constructed in Theorem \ref{th}
 does not belong to the set
 $$\mathrm{Bad} = \{(\theta_1,\theta_2) | \inf_{x \in \mathbb{N}} x^{1/2} \max (||\theta_1x||,||\theta_2x||)>0\}.$$

{\bf Remark 2.}
Let $\pmb{\theta} = \left(\dfrac{a_{1}}{q}, \dfrac{a_{2}}{q}\right)$ be rational.
Let $\pmb{\eta} = (\eta_1,\eta_2) \notin \dfrac{1}{q} \cdot \mathbb{Z}^2$, then for any $x \in \mathbb{Z}$,
$$
\max_{i=1,2} \left|\left|x \dfrac{a_{i}}{q}-\eta_i\right|\right| \ge {\rm dist} \left(\pmb{\eta},\dfrac{1}{q} \cdot \mathbb{Z}^2\right)> 0 .
$$
So the set $$\mathcal{B} = \left\{ \pmb{\eta} :  \inf_{x \in \mathbb{Z}} \max_{i=1,2} \left|\left|x \dfrac{a_{i}}{q}-\eta_i\right|\right| > 0 \right\} $$
contains $\mathbb{R}^2 \setminus \dfrac{1}{q} \cdot \mathbb{Z}^2 $ and is trivially winning.
It is clear, that for any one-dimensional affine subspace $\ell$ we have $\mathcal{B} \cap \ell \supset \left(\mathbb{R}^2 \setminus \dfrac{1}{q} \cdot \mathbb{Z}^2 \right)\cap \ell$.
So obviously $\mathcal{B} \cap \ell$ is also winning in $\ell.$
\section{linear dependent case}\label{dep}
Let $1,\theta_1,\theta_2$ be linear dependent and at least one of $\theta_j$ is irrational.  It means that there exists $\bm{z} = (z_0,z_1,z_2) \in \mathbb{Z}^3$ such that $(\bm{z} , \bm{\theta}) = 0.$
Let us consider two-dimensional rational subspace $$\pi = \{\pmb{x} \in R^3 : (\pmb{x},\pmb{z}) =0\},$$ so $\bm{\theta} \in \pi$.

Let us define one-dimensional subspace $\mathcal{P} = \{ (x_1, x_2) : (1,x_1,x_2) \in \pi\} \subset \mathbb{R}^2$ .

We will prove that there exists a constant $\gamma$ such that for any
$
\eta = (\eta_1,\eta_2) \in \mathcal{P}
$
the inequality
$$
\qquad \max_{i=1,2} ||\theta_i x - \eta_i|| < \frac{\gamma}{x}
$$
has  infinitely many solutions in $x \in \mathbb{N}$.
(This statement is similar to Chebushev's theorem \cite[Theorem 24, Chapter 2]{CH}

Denote by $\Lambda = \pi \cap \mathbb{Z}^3$ the integer lattice with the determinant $ d : = \det \Lambda = |\pmb{z}|.$
Denote by $\{\pmb{g}_{\nu} = (q_{\nu}, a_{1 \nu}, a_{2 \nu})\}_{\nu=1,2,3 \ldots} \subset \Lambda$ the sequence of the best approximations of $\pmb{\theta}$ by the lattice $\Lambda$ and corresponding parallelograms

$$
\Pi_{\nu} = \{ \pmb{x} = (x_0,x_1,x_2) \in \pi : 0\le x_0 \le q_{\nu} : $$
$${\rm dist} (\pmb{x}, l(\pmb{\theta})) \le {\rm dist} (\pmb{g}_{\nu -1}, l(\pmb{\theta})) \}
$$
which contains a fundamental domain of the two-dimensional $\Lambda.$
Obviously, vol $\Pi_{\nu} \le 4d$.
So
\begin{equation}\label{1}
 {\rm dist} (\pmb{g}_{\nu -1}, l(\theta)) \ll \dfrac{d}{q_{\nu}}
\end{equation}
with an absolute constant in the sign $\ll$.
 It is clear that for any point $\bm{\eta} \in \pi$ the shift  $  \bm{\eta} + \Pi_{\nu}$ contains a point of $\Lambda$.

For any $\bm{\eta} = (\eta_1,\eta_2) \in \mathcal{P}$ and for any positive integer $\nu$ the planar domain $\overline{\pmb{\eta}}+ \Pi_{\nu}$, $\overline{\pmb{\eta}} = (1, -\eta_1, -\eta_2)$ contains an integer point
$\bm{y} = (x,y_1,y_2) \in \Lambda$.

It is clear that
\begin{equation}\label{2}
1\le x \le 1 + q_{\nu}
\end{equation}
and
$$
\max_{i=1,2}
 || \theta_i x- \eta_i||  \ll {\rm dist} (\pmb{y},l(\pmb{\theta}) + \overline{\pmb{\eta}}) \ll  {\rm dist} (l(\pmb{\theta}),\pmb{g}_{\nu -1}),
$$
and by \eqref{1}
\begin{equation}\label{3}
\max_{i=1,2} || \theta_i x- \eta_i|| \ll \dfrac{d}{q_{\nu}}
\end{equation}

From \eqref{2}, \eqref{3} it follows that the inequality
$$ \max_{i=1,2} || \theta_i x- \eta_i|| \ll \dfrac{d}{x}.$$
has infinitely many solutions
and everything is proved.

\section{Inductive construction of integer points}\label{indep}
Let $\omega(t)$ be arbitrary positive function monotonically (slowly) increasing to infinity. Here we describe the inductive construction of integer points
$\pmb{z}_{\nu}=(q_{\nu},z_{1 \nu },z_{2 \nu })$.
The base of the induction process is trivial. One can take an arbitrary primitive pair of integer vectors that can be completed to a basis of $\mathbb{Z}^3$.

Suppose that we have two primitive integer vectors
$$
\pmb{z}_{\nu-1}=(q_{\nu-1},z_{1 \, \nu-1 },z_{2 \, \nu-1 })\in \mathbb{Z}^3
$$
$$
\pmb{z}_{\nu}=(q_{\nu},z_{1 \, \nu},z_{2 \, \nu})\in \mathbb{Z}^3.
$$
Now we explain how to construct the next integer vector $\pmb{z}_{\nu+1}.$

We consider two-dimensional subspace
$$\pi_{\nu} = <\pmb{z}_{\nu-1},\pmb{z}_{\nu}>_{\mathbb{R}}. $$
As the pair of vectors $\pmb{z}_{\nu-1}$ and $\pmb{z}_{\nu}$  is primitive, the lattice
$$\Lambda_{\nu} := <\pmb{z}_{\nu-1},\pmb{z}_{\nu}>_{\mathbb{Z}}=  \pi_{\nu}\cap \mathbb{Z}^3 .$$
By $d_{\nu} = \det \Lambda_{\nu}$ we denote the two-dimensional fundamental volume of the lattice $\Lambda_{\nu}.$
Now we define vector $\pmb{n_{\nu}} = (n_{0 \nu},n_{1\nu}, n_{2\nu})\in \mathbb{R}^3$ from the conditions
$$\pi_{\nu}= \{\pmb{x} \in \mathbb{R}^3 : (\pmb{x},\pmb{n_{\nu}}) =0 \}, \quad |\pmb{n}_{\nu}|=1.$$
Put
\begin{equation}\label{sigma}
 \sigma_{\nu}= {\rm dist} (\pmb{z}_{\nu -1}, l(\pmb{z}_{\nu}))
\end{equation}
Obviously, $|\pmb{z}_{\nu}| \asymp q_{\nu}$ and
\begin{equation}\label{asymp}
\sigma_{\nu}\asymp \dfrac{d_{\nu}}{q_{\nu}}.
\end{equation}
We define a vector $\pmb{e}_{\nu}$
from the conditions
\begin{equation}\label{e}
\pmb{e}_{\nu} \in \pi_{\nu} , \qquad |\pmb{e}_{\nu}| = 1, \qquad (\pmb{e}_{\nu}, \pmb{z}_{\nu}) = 0
\end{equation}
so $\pmb{e}_{\nu}$ is parallel to $\pi_{\nu}$ and orthogonal to $\pmb{z}_{\nu}$.

Define the rectangle $$\Pi_{\nu}
= \{ \pmb{x} = (x_0,x_1,x_2)  : \pmb{x} = t\pmb{z}_{\nu} + r\pmb{e}_{\nu}, \,\, 0 \le t\le |\pmb{z}_{\nu}|, \,\, |r| \le \sigma_{\nu}  \}.$$
It is clear that rectangle $\Pi_{\nu} \subset \pi_{\nu}$ contains a fundamental domain of the lattice $\Lambda_{\nu}.$
We need two axillary vectors $\pmb{z}^a_{\nu}$ and $\pmb{z}^b_{\nu}$ defined as
$$\pmb{z}^a_{\nu} = \pmb{z}_{\nu} + a_{\nu} \pmb{e}_{\nu}$$
$$
\pmb{z}^{b}_{\nu} = \pmb{z}^a_{\nu} + b_{\nu} \pmb{n}_{\nu}
$$
where positive $a_{\nu}$ is chosen in such a way
\begin{equation}\label{a}
a_{\nu} d_{\nu}^2 \le \nu^{-1} \omega \left( \dfrac{q_{\nu}^2}{d_{\nu}^2} \cdot \dfrac{1}{a_{\nu}} \right)
\end{equation}
and
\begin{equation}\label{b}
b_{\nu} = a_{\nu} {\rm min} \left(1, \dfrac{d_{\nu}}{q_{\nu}}\right).
\end{equation}

From the construction it follows that
\begin{equation}\label{|Z|}
|\pmb{z}_{\nu}^a| \asymp |\pmb{z}_{\nu}^b| \asymp |\pmb{z}_{\nu}| \asymp q_{\nu}.
\end{equation}
Integer lattice $\mathbb{Z}^3$ splits into levels with respect to the two-dimensional sublattice $\Lambda_{\nu}$ in such a way that
$$\mathbb{Z}^3 = \bigsqcup_{i\in \mathbb{Z}} \Lambda_{\nu,i},$$
where $\Lambda_{\nu,j}=\Lambda_{\nu} + j \pmb{z'},  j\in \mathbb{Z}$ and integer vector $\pmb{z'}$ completes the couple $\pmb{z}_{\nu-1}, \pmb{z}_{\nu}$ to the basis in $\mathbb{Z}^3$.
We consider affine subspace
$\pi^1_{\nu}= \pi_{\nu} + \pmb{z}' \supset \Lambda_{\nu,1},$ which is parallel to $\pi_{\nu}$.
It is clear that ${\rm dist} (\pi_{\nu},\pi^1_{\nu}) = \frac{1}{d_{\nu}}$.

We need to determine the next integer point $z_{\nu+1}$.
Denote by $\mathfrak{P}$ central projection with center 0 onto the affine subspace
 $\pi^1_{\nu}$.
  We consider triangle $\Delta$ with the vertices $z_{\nu},z^a_{\nu},z^{b}_{\nu}$ and its image $\mathfrak{P} \Delta$ under the projection $\mathfrak{P}$.
   Define \begin{equation}\label{Z}
  \pmb{Z} = \mathfrak{P} {\pmb{z}^b_{\nu}}.
   \end{equation}
    One can see that
  \begin{equation}\label{|zz|}
  |\pmb{Z}| \asymp \dfrac{q_{\nu}}{d_{\nu} b_{\nu}}.
  \end{equation}
  \includegraphics[width=10cm,height=5cm]{fig1.jpg}
  \captionof{figure}{The central projection $ \mathfrak{P}$.}\label{fig:intab}

Define rays
$$ \mathcal{R}_1 = \{\pmb{z} = \pmb{Z} + t \pmb{z}_{\nu}, \quad t\ge 0\} \quad \text{and} \quad \mathcal{R}_2 = \{\pmb{z} = \pmb{Z} + t \pmb{z}^a_{\nu}, \quad t\ge 0\}.$$
It is clear that $\mathcal{R}_1 \cap \mathcal{R}_2 = \{\pmb{Z}\}$ and $\mathcal{R}_1, \mathcal{R}_2 \subset \pi^1_{\nu}$.
Moreover, the whole convex angle bounded by rays $\mathcal{R}_1, \mathcal{R}_2$ form the image of the triangle $\Delta$ under the projection $\mathfrak{P}$:
$$
\mathfrak{P}\Delta = {\rm conv}\, (\mathcal{R}_1 \cup \mathcal{R}_2) .
$$
Affine subspace $\pi_{\nu}^1$ contains the affine lattice $\Lambda_{\nu}^1 = \Lambda_{\nu} + \pmb{z}'$ which is congruent to the lattice $\Lambda_{\nu}.$
Thus, for any $\pmb{\zeta} \in \pi_{\nu}^1$ the shift $\Pi_{\nu}+\pmb{\zeta}$ contains an integer point from $\Lambda_{\nu}^1.$

Put
\begin{equation}\label{tau}
\tau_{\nu} = \dfrac{2\sigma_{\nu} |\pmb{z}_{\nu}|}{a_{\nu}},
\end{equation}
consider the point
$$\pmb{\zeta}_{\nu} =  \pmb{Z}+ \tau_{\nu} \pmb{z}_{\nu} +\sigma_{\nu}\pmb{e}_{\nu}\in \pi_{\nu}^1. $$
and the rectangle
$$
\Pi_{\nu}^1 = \Pi_{\nu} +\pmb{\zeta}_{\nu} \subset \pi_{\nu}^1.
$$
It is clear that
$$
\Pi_{\nu}^1 \subset \mathfrak{P} \Delta$$
(here $\pmb{Z}$ was defined in \eqref{Z}, $\pmb{e}_{\nu}$ was defined in \eqref{e} and parameters $\sigma_{\nu}, \tau_{\nu}$ come from \eqref{sigma} and \eqref{tau}).

Now we take the integer point
$$\bm{z}_{\nu+1} = (q_{\nu+1},z_{1 \, \nu+1 },z_{2 \, \nu+1 }) \in \Lambda_{\nu}^1 \cap \Pi_{\nu}^1.$$
From the construction it follows that
$$q_{\nu+1} \asymp |\pmb{z}_{\nu+1}| \asymp  |\pmb{Z}|+ \tau_{\nu}  |\pmb{z}_{\nu}| + |\pmb{z}_{\nu}|
$$
$$
\asymp q_{\nu} \left( 1+ \dfrac{1}{d_{\nu} b_{\nu}} + \dfrac{\sigma_{\nu}}{a_{\nu}} \right)  \asymp
 q_{\nu} \left( 1+ \dfrac{1}{d_{\nu} b_{\nu}} \right)+ \dfrac{d_{\nu}}{a_{\nu}} \asymp \dfrac{q_{\nu}}{d_{\nu} b_{\nu}}.
$$(Here we use (\ref{asymp},\ref{b}, \ref{|Z|}, \ref{|zz|},  \ref{tau}).)
From \eqref{b} we see that
\begin{equation}\label{q+1}
q_{\nu+1} \gg \left( \dfrac{q_{\nu}}{d_{\nu}}\right)^2 \dfrac{1}{a_{\nu}}.
\end{equation}
Now we are able to define the next  two-dimensional lattice
$$
\Lambda_{\nu+1} = <\pmb{z}_{\nu},\pmb{z}_{\nu+1}>_{\mathbb{Z}}.
$$
Let $d_{\nu+1}$ be its fundamental volume. We will estimate the value of $d_{\nu+1}$ taking into account \eqref{b}
\begin{equation}\label{d}
d_{\nu+1} \ll q_{\nu} \cdot {\rm dist} (\pmb{z}_{\nu+1}, l(\pmb{z}_{\nu})) \ll \dfrac{q_{\nu}}{d_{\nu}} \cdot \dfrac{a_{\nu}}{b_{\nu}} \ll \left(\dfrac{q_{\nu}}{d_{\nu}} \right)^2 \ll q_{\nu}^2.
\end{equation}
From \eqref{q+1} and \eqref{d} we deduce that
$$
d_{\nu+1}  \ll a_{\nu}  d_{\nu}^2 q_{\nu+1}.
$$
By the choice of $a_{\nu}$ (by formula \eqref{a}) we have
\begin{equation}\label{7}
d_{\nu+1} \le \dfrac{\omega(q_{\nu+1})}{\nu}.
\end{equation}
\section{The vector $\pmb{\theta}$}

Now we define
$$
\pmb{\theta}_{\nu} = (\theta_{1 \nu},\theta_{2 \nu}), \quad \theta_{j \nu} = \dfrac{q_{j \nu}}{q_{\nu}}.
$$
We consider the angles between the successive vectors $\pmb{n}_{\nu}$ and $\pmb{n}_{\nu+1}$:
$$
\alpha_{\nu} ={\rm angle} (\pmb{n}_{\nu}, \pmb{n}_{\nu+1}) \asymp \tan {\rm angle} (\pmb{n}_{\nu}, \pmb{n}_{\nu+1})
$$
 \includegraphics[width=15cm,height=5cm]{pic21.jpg}
 \captionof{figure}{The vector $\pmb{z}_{\nu+1}$ intersects the interior of the triangle
$\Delta = \pmb{z}_{\nu} \pmb{z}_{\nu}^a\pmb{z}_{\nu}^b$.}\label{fig:intab}

As $\pmb{z}_{\nu+1} \in \mathfrak{P} \Delta$ (see Fig. 2), we have
$$
\tan  {\rm angle} (\pmb{n}_{\nu}, \pmb{n}_{\nu+1}) \le \dfrac{b_{\nu}}{a_{\nu}},
$$
and so
\begin{equation}\label{alpha}
\alpha_{\nu} \ll \dfrac{b_{\nu}}{a_{\nu}}.
\end{equation}
As $\pmb{z}_{\nu+1} \in \mathfrak{P} \Delta$ we have
\begin{equation}\label{theta}
|\pmb{\theta}_{\nu} - \pmb{\theta}_{\nu+1}| \ll \dfrac{\sqrt{a_{\nu}^2 + b_{\nu}^2}}{q_{\nu}} \ll \dfrac{a_{\nu}}{q_{\nu}},
\end{equation}
by the same argument.
There exist limits $$\lim_{\nu \rightarrow \infty} \pmb{\theta}_{\nu} = \pmb{\theta} = (\theta_1,\theta_2) \quad {\rm and} \quad \lim_{\nu \rightarrow \infty} \pmb{n}_{\nu} = \pmb{n},$$
and from (\ref{alpha}) and (\ref{theta}) we deduce that
\begin{equation}\label{5}
0 < |\pmb{\theta} - \pmb{\theta}_{\nu}| \ll \dfrac{a_{\nu}}{q_{\nu}}
\end{equation}
and
\begin{equation}\label{angle}
{\rm angle} (\pmb{n},\pmb{n}_{\nu}) \ll \dfrac{b_{\nu}}{a_{\nu}}
\end{equation}
It is clear that $\pmb{\theta} \not\in \mathbb{Q}^2.$
A slight modification\footnote{A similar procedure was explained in \cite{NGM1}. There author provides the linear independence of coordinates of the limit vector by "going away from all rational subspaces" (the beginning of proof of Theorem 1 in the case $k=1$, p.132 and the beginning of section 5, p.146).} of the procedure of choosing vectors $\pmb{z}_{\nu}$ ensures the condition that $1,\theta_1,\theta_2$ are linear independent over $\mathbb{Z}$. Define $\pi = \{\pmb{x} \in \mathbb{R}^3: (\pmb{x},\pmb{n}) =0\}$. Then $\pmb{\theta} \in \pi$ by continuity and we can assume that $\pmb{n} \not\in \mathbb{Q}^3.$
\section{Winning property}
Consider  one-dimensional affine subspaces  $$\mathcal{P}_{\nu} = \{(x_1,x_2) \in \mathbb{R}^2 : (1,x_1,x_2) \in \pi_{\nu}\} \subset \mathbb{R}^2$$ and  $$\mathcal{P} = \{(x_1,x_2) \in \mathbb{R}^2 : (1,x_1,x_2) \in \pi\} \subset \mathbb{R}^2,$$ where $\pi$ was defined in the end of the previous section.
Let
$$ B_1(\pmb{\theta}) = \{ \pmb{\xi} \in \mathbb{R}^2 : {\rm dist} (\pmb{\xi},\pmb{\theta}) < 1\}. $$
We will show that for any $\pmb{\eta} = (\eta_1, \eta_2) \in \mathcal{P} \cap B_1 (\pmb{\theta})$ there exists infinitely many solutions of the inequality
$$
\max_{i=1,2}
 || \theta_{i} x- \eta_{i}|| < \dfrac{\omega(x)}{x}
$$
in integer $x.$
Denote by $\pmb{\eta}_{\nu} = (\eta_{1 \nu},\eta_{2 \nu})$ the orthogonal projection of $\pmb{\eta}$ onto $\mathcal{P}_{\nu}$.
From \eqref{angle} we see that
\begin{equation}\label{4}
|\pmb{\eta} -\pmb{\eta}_{\nu}| \ll \dfrac{b_{\nu}}{a_{\nu}}
\end{equation}

For any $\pmb{\eta}_{\nu} = (\eta_{1 \nu}, \eta_{2 \nu}) \in \mathcal{P}_{\nu}$
 the planar domain $\overline{\pmb{\eta}_{\nu}}+ \Pi_{\nu}, \overline{\pmb{\eta}_{\nu}} = (1, -\eta_{1 \nu}, -\eta_{2 \nu})$ contains an integer point
$\pmb{y}_{\nu} = (x_{\nu},y_{1 \nu},y_{2 \nu}) \in \Lambda_{\nu}$.
It is clear that
\begin{equation}\label{12}
|x_{\nu}|\ll q_{\nu}
\end{equation}
and
\begin{equation}\label{13}
\max_{i=1,2}
 | \theta_{i \nu} x_{\nu}- \eta_{i \nu}-y_{i \nu}|  \ll \dfrac{d_{\nu}}{q_{\nu}}.
\end{equation}
By \eqref{5},\eqref{4},\eqref{12},\eqref{13} we have
\begin{equation*}
\max_{i=1,2}
 || \theta_{i} x_{\nu}- \eta_{i}|| \le
  |x_{\nu}| \max_{i=1,2}  |\theta_i - \theta_{i \nu}| +  \max_{i=1,2} ||\theta_{i \nu} x_{\nu}- \eta_{i \nu}|| + \max_{i=1,2} |\eta_i - \eta_{i \nu}|
\end{equation*}
 $$
 \ll a_{\nu} + \dfrac{d_{\nu}}{q_{\nu}} + \dfrac{b_{\nu}}{a_{\nu}} \ll \dfrac{d_{\nu}}{q_{\nu}}
$$
in the last inequality we use \eqref{b}.
By \eqref{7} we have
$$
\max_{i=1,2}
 || \theta_{i} x_{\nu}- \eta_{i}|| \le \dfrac{\omega(q_{\nu})}{q_{\nu}}
$$
for large $\nu.$
As $\pmb{\overline{\eta}} \in \pi$ and $\pmb{y}_{\nu} \in \pi_{\nu}$,
$\max_{i=1,2}
 || \theta_{i} x_{\nu}- \eta_{i}|| \neq 0$ infinitely often (in fact for all large $\nu$.)