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

About 10 per cent of A and B main sequence stars possess stable, globally organised magnetic fields with a strength of at least 300~G \citep{auriere:2007,sikora:2019a}. These stars typically exhibit anomalous absorption spectra, shaped by surface overabundances of Si, Fe-peak, and rare-earth elements, and are known as magnetic chemically peculiar or ApBp stars. The process of radiatively driven atomic diffusion \citep{michaud:2015} responsible for these abundance anomalies also produces a high-contrast, long-lived non-uniform horizontal distribution of chemical elements. Local variation of metal abundances associated with these chemical spots modifies emergent stellar radiation and leads to a characteristic photometric rotational modulation known as $\alpha^2$~CVn (ACV) type of stellar variability \citep{samus:2017}.

The properties of magnetic fields of ApBp stars do not depend on stellar mass or rotation rate, making them strikingly different from characteristics of the dynamo-generated fields observed in late-type stars \citep[e.g.][]{vidotto:2014}. It is believed that the strong fields found in early-type stars are dynamically stable, `fossil' remnants of the magnetic flux generated or acquired by these stars at some earlier evolutionary phase \citep{braithwaite:2004,neiner:2015}. The origin of fossil fields is not understood. This magnetic flux might be inherited from molecular clouds at stellar birth \citep{mestel:1999}, produced by the convective dynamo during pre-main sequence evolution \citep{moss:2004} or created by a short-lived dynamo operating during stellar merger events \citep{schneider:2019}.

The binary characteristics of early-type magnetic stars may provide crucial clues, allowing one to test alternative fossil field hypotheses. The non-magnetic chemically peculiar stars of Am (A-type stars with enhanced lines of Fe-peak elements) and HgMn (late-B stars identified by strong lines of Hg and/or Mn) types are frequently found in close binaries (\citealt*{gerbaldi:1985}; \citealt{ryabchikova:1998a}; \citealt{carquillat:2007}), including eclipsing systems \citep{nordstrom:1994,strassmeier:2017,takeda:2019}. In contrast, only about ten close ($P_{\rm orb}<20$~d) spectroscopic binaries containing at least one magnetic ApBp star are known \citep{landstreet:2017}. The overall incidence rate of magnetic upper main sequence stars in close binaries is less than 2 per cent \citep{alecian:2015a}, although this fraction is significantly higher if one includes wide long-period systems \citep{mathys:2017}. This low incidence of magnetic ApBp stars in close binaries is frequently considered as an argument in favour of the stellar merger origin of fossil fields \citep{de-mink:2014,schneider:2016}. In this context, confirmation of magnetic ApBp stars in short-period binary systems gives support to alternative theories or, at least, demonstrates that early-type stars may acquire magnetic fields through different channels. In addition, detached close binary stars, particularly those showing eclipses, are valuable astrophysical laboratories that provide model independent stellar parameters and allow one to study pairs of co-evolving stars formed in the same environment. Until recently, no early-type magnetic stars in eclipsing binaries were known. The first such system, HD\,66051, was identified by \citet{kochukhov:2018b}. The second system, HD\,62658 containing twin components of which only one is magnetic, was found by \citet{shultz:2019}. Several other eclipsing binaries containing candidate ApBp stars were proposed (\citealt{hensberge:2007}; \citealt*{gonzalez:2010a}; \citealt{skarka:2019}), but the magnetic nature of these stars has not been verified by direct detections of their fields using the Zeeman effect. In this paper we put a spotlight on another candidate eclipsing magnetic Bp star, which received little attention prior to our work despite being significantly brighter than the confirmed magnetic eclipsing systems HD\,62658 and HD\,66051.

V772~Cas\ (HR\,481, HD\,10260, HIP\,7939) is a bright ($V=6.7$) but little studied chemically peculiar late-B star. The exact type of its spectral peculiarity is uncertain. \citet{cowley:1972} classified this star as B8IIIpSi whereas \citet{dworetsky:1976} considered it an HgMn star. Both studies used low-dispersion classification spectra. The former BpSi classification appears to be more common in recent literature \citep[e.g.][]{gandet:2008,renson:2009,skarka:2019}. The variable star designation comes from \citet{kazarovets:1999}, who suggested this object to be an $\alpha^2$~CVn-type variable based on the Hipparcos epoch photometry \citep{perryman:1997}. \citet{otero:2007} discovered eclipses in the Hipparcos light curve. This analysis was improved by \citet{gandet:2008}. He confirmed the presence of primary eclipses, derived an orbital period of 5.0138~d and demonstrated that archival photographic radial velocity measurements \citep{hube:1970} show coherent variation with the same period. These results, along with the B8IIIpSi spectral classification and an evidence of the out-of-eclipse variability, led \citet{gandet:2008} to suggest V772~Cas\ as an $\alpha^2$~CVn variable in a short-period eclipsing binary -- an exceptionally rare and interesting object akin to the recently discovered magnetic eclipsing binaries HD\,66051 and HD\,62658. Apart from the study by \citet*{huang:2010}, who determined $T_{\rm eff}$\,=\,$13188\pm250$~K and $\log{g}$\,=\,$3.43\pm0.05$ from the low-resolution hydrogen Balmer line spectra, no model atmosphere and/or abundance analysis was carried out for V772~Cas. In fact, to the best of our knowledge, this star was never studied with high-resolution spectra.

In this paper we present a detailed photometric and spectroscopic investigation of V772~Cas\ that provides a new insight into the nature of this star. In Sect.~\ref{sec:obs} we describe the new observational data used in our study. Orbital modelling and derivation of the binary component parameters are presented in Sect.~\ref{sec:binary}. Model atmosphere parameters and chemical abundances of the primary are determined in Sect.~\ref{sec:primary}. The paper concludes with the discussion in Sect.~\ref{sec:disc}.

\section{Observations}
\label{sec:obs}

\subsection{Space photometry}

The NASA Transiting Exoplanet Survey Satellite \citep[TESS;][]{ricker:2015} began its nominal 2 year mission in 2018 to discover Earth-sized transiting exoplanets. Using four cameras that cover a combined field of view of 24$^{\circ}$ $\times$ 96$^{\circ}$ and a red filter that records light in the range 6000 to 10500\AA, both ecliptic hemispheres are surveyed for one year each, in 13 sectors that extend from the ecliptic plane to the pole. Each sector is observed for approximately 27.5 days, and a given field on the sky can be observed in multiple sectors if it falls on overlap regions. Nominal TESS targets are bright, with $I_c$ magnitudes between 4 and 13, with a noise floor of approximately 60 ppm hr$^{-1}$.

Full Frame Images (FFIs) from TESS are available at a 30-minute cadence for the entire field of view, allowing light curves to be extracted for all objects that fall on the detector. Certain high priority targets were pre-selected by the TESS mission to be observed with 2-minute cadence, some of which were chosen from guest investigator programs.

V772 Cas was observed in cycle 2 of the TESS mission in sector 18 (2019-Nov-02 to 2019-Nov-27) in 30-minute cadence mode (and was not pre-selected for 2-minute cadence observations; TESS Input Catalog ID: 444833007). A target pixel file of a 50 $\times$ 50 grid of pixels centered on V772 Cas was downloaded with TESScut\footnote{\url{https://mast.stsci.edu/tesscut/}} \citep{Brasseur2019}, with further processing aided by the \textsc{lightkurve} package \citep{Lightkurve2018}.
A light curve for V772~Cas\ was extracted with aperture photometry using an aperture of 36 pixels centred on the target. This aperture was chosen to minimise contaminating flux from neighbouring stars, which would otherwise bias the binary system modelling and possibly introduce additional signals from other sources, while still achieving a high signal-to-noise ratio. A principal component analysis routine was then applied to the extracted light curve to detrend against signals common to neighbouring pixels (outside of a 15 $\times$ 15 pixel exclusion zone). After removing outliers and data points more prone to systematic effects (mostly scattered light near TESS orbital perigee), the light curve includes 958 observations spanning 22.1 days.

The Hipparcos light curve analysed by \citet{otero:2007} and \citet{gandet:2008} covered 233 orbital cycles with about a dozen photometric data points tracing the primary eclipse. Although the TESS data corresponds to a time span of just 4.4 orbits, it samples the eclipses much more densely owing to its higher measurement cadence. Moreover, individual TESS measurements are about a factor of 200 more precise than the Hipparcos photometry.

\subsection{High-resolution spectroscopy}

We obtained high-resolution spectra of V772~Cas\ using the High-Efficiency and high-Resolution Mercator \'Echelle Spectrograph (HERMES) mounted on the 1.2 m Mercator telescope at the Observatorio del Roque de los Muchachos, La Palma, Canary Islands, Spain. This instrument provides coverage of the 3700--9100~\AA\ wavelength region at a resolving power of 85\,000 \citep{raskin:2011}. The star was observed on five consecutive nights, from January 27 to January 31, 2020, with two spectra obtained on each night. Exposure times of 600--800~s were used, yielding a signal-to-noise ratio ($S/N$) of 190--280 per $\approx$\,0.03~\AA\ pixel of the extracted spectrum in the wavelength interval from 5000 to 5500~\AA. The HERMES pipeline reduction software was employed to perform the basic \'echelle data reduction steps, including the bias and flat field corrections, extraction of one-dimensional spectra and wavelength calibration. The resulting merged, un-normalised spectra were then normalised to the continuum with the methodology described by \citet{rosen:2018}.

The log of ten HERMES observations of V772~Cas\ is given in Table~\ref{tab:obs}. The first four columns of this table list the UT date of observation, the corresponding heliocentric Julian date, the orbital phase calculated according to the ephemeris derived in Sect.~\ref{sec:binary}, and the $S/N$ ratio. The last column provides radial velocities determined in the next section.

\subsection{Spectroscopic classification and radial velocity measurements}

Initial qualitative analysis of the high-resolution spectra of V772~Cas\ showed the presence of a single set of spectral lines with variable radial velocity and line strengths typical of HgMn late-B chemically peculiar stars. Specifically, the HgMn classification of this star is unambiguously demonstrated by the presence of \ion{Hg}{ii} 3984~\AA, \ion{Ga}{ii} 6334~\AA\ and numerous strong \ion{Mn}{ii} and \ion{P}{ii} lines. On the other hand, ionised Si lines are not anomalously strong and rare-earth absorption features are not prominent in the spectrum of V772~Cas, which argues against identification of this object as a Si-peculiar Bp star.

Based on this initial assessment, we extracted a line list from the VALD database (\citealt{ryabchikova:2015}; \citealt*{pakhomov:2019}) with the chemical abundances typical of HgMn stars \citep*{ghazaryan:2018} and model atmosphere parameters $T_{\rm eff}$\,=\,13000~K, $\log{g}$\,=\,3.5, close to the values reported by \citet{huang:2010}\footnote{Precise choice of atmospheric parameters and abundances is not important for the multi-line method used in this paper.}. This line list was employed for calculation of least-squares deconvolved (LSD) profiles \citep*{kochukhov:2010a} with the goal to study the radial velocity variation of the primary and search for spectral signatures of the secondary. The LSD profiles were constructed by combining 1094 metal lines deeper than 5 per cent of the continuum. These average spectra are illustrated in Fig.~\ref{fig:lsd}, where the data are phased with the orbital ephemeris $HJD=2458803.4016+5.0138\times E$ determined below (Sect.~\ref{sec:binary}).

The radial velocity of the primary was measured from the LSD profiles with the help of the centre-of-gravity method. These measurements made use of the $\pm30$~km\,s$^{-1}$\ velocity range around the line centre.  The resulting radial velocity changes from $-37$ to 38~km\,s$^{-1}$. This variation is coherent on the 5~d time-scale covered by the spectroscopic observations and is consistent with the orbital period seen in photometry. Individual radial velocities are reported in the last column of Table~\ref{tab:obs}. The formal uncertainty of these measurements is 50--60~m\,s$^{-1}$.

At the same time, no evidence of a spectral contribution of the secondary was found in the LSD profiles. This is not surprising considering the luminosity ratio of $>$\,100 estimated from the eclipse depths in the TESS light curve. We tested alternative solar-composition LSD line masks for $T_{\rm eff}$\ of 7000, 9000, and 11000~K but still could not detect any features in the LSD profile time series that could be attributed to the secondary component.

\begin{table}
\caption{Log of spectroscopic observations of V772~Cas. \label{tab:obs}}
 \begin{tabular}{llllr}
\hline\hline
UT date & HJD & Orbital & $S/N$ & $V_{\rm r}~~~~$ \\
        &     & phase & (pixel$^{-1}$)           & (km\,s$^{-1}$) \\
\hline
2020-01-27 & 2458876.3638 & 0.552 & 256 &  $ 11.17\pm0.05$ \\
2020-01-28 & 2458876.5045 & 0.580 & 237 &  $ 17.54\pm0.06$ \\
2020-01-28 & 2458877.3365 & 0.746 & 223 &  $ 37.70\pm0.06$ \\
2020-01-29 & 2458877.5078 & 0.780 & 256 &  $ 36.87\pm0.06$ \\
2020-01-29 & 2458878.3251 & 0.943 & 237 &  $ 11.91\pm0.06$ \\
2020-01-30 & 2458878.5005 & 0.978 & 238 &  $  3.91\pm0.06$ \\
2020-01-30 & 2458879.3307 & 0.144 & 277 &  $-32.46\pm0.05$ \\
2020-01-31 & 2458879.5036 & 0.179 & 255 &  $-37.10\pm0.06$ \\
2020-01-31 & 2458880.3289 & 0.343 & 232 &  $-34.51\pm0.06$ \\
2020-02-01 & 2458880.5186 & 0.381 & 191 &  $-28.49\pm0.06$ \\
\hline
\end{tabular}
\end{table}

\begin{figure}
\centering
\figps{0.85\hsize}{0}{lsd_profiles}
\caption{Least-squares deconvolved profiles of V772 Cas. Profiles for different orbital phases are shifted vertically with a step of 0.2. The dotted lines show continuum level for each spectrum. The solid lines show LSD profile for the entire velocity span with the thicker segments indicating the velocity internal employed for radial velocity measurements. The orbital phase corresponding to each observation is given on the right.}
\label{fig:lsd}
\end{figure}

\section{Binary system modelling}
\label{sec:binary}

We subject the TESS photometry and radial velocity measurements to simultaneous modelling to
determine the fundamental parameters of the components of V772 Cas. To carry out this simultaneous
modelling, we wrap the {\sc PHOEBE} binary modelling code \citep{prsa:2005,prsa:2011} into the Markov
Chain Monte Carlo (MCMC) ensemble sampling code {\sc emcee} \citep{foreman-mackey:2013}. This
methodology has been outlined and employed in the modelling of the magnetic systems HD~66051
\citep{kochukhov:2018b} and HD~62658 \citep{shultz:2019} and is briefly described below.

As discussed in detail by \citet{prsa:2005} and \citet{prsa:2016}, the {\sc PHOEBE} code employs a generalised Roche geometry model of binary systems which treats the orbital motion of eclipsing binary-star components along with the surface brightness variation caused by limb darkening, gravity darkening, reflection, ellipsoidal variations due to non-spherical shapes of the components, and brightness spots on their surfaces. The two latter phenomena are of particular interest here as possible explanations of the out-of-eclipse modulation observed in the TESS light curve of V772~Cas. The ellipsoidal variability arises due to the changing cross-section size that faces the observer. It occurs with the orbital period and has a distinctive shape, with two maxima and minima per orbital cycle. The minimal light corresponds to eclipse phases when the components are aligned with the line of sight. In contrast, the flux variation associated with surface spots occurs with rotational periods of the components, which are not necessarily the same as the orbital period. The shapes and amplitudes of spot-induced light curves are diverse, with the brightness extrema generally not linked to particular orbital phases. In the context of {\sc PHOEBE} analysis ellipsoidal variation is an integral part of binary system modelling whereas spots are introduced with a set of additional free parameters.

\subsection{Modelling setup and results}
\label{sec:binary_results}

Given a set of input parameters {\sc PHOEBE} produces a forward binary model that consists
of both a photometric model and radial velocity model. The aim of our modelling is to
obtain the input parameters which produce the {\sc PHOEBE} model that best reproduces the
entire observed TESS light curve, including eclipses and out-of-eclipse variation, and HERMES radial velocities. After we determine a
reasonable starting model by hand, we use {\sc emcee} to explore the posterior
distributions of the free model parameters through ensemble MCMC sampling. We use 128
chains and initially run the code until it has reached 1000 iterations beyond convergence
defined by a less than 1 per cent change in the estimated autocorrelation time \citep{foreman-mackey:2019}.
After convergence is reached, we record the final positions of each chain, discard all
previous steps as burn-in, and re-initialize the algorithm from these points for an
additional 5000 iterations. This results in 640\,000 model evaluations. We marginalize
over all varied parameters to arrive at posterior distributions for these parameters. From
these we draw the median and 68 per cent highest posterior density confidence intervals as the
parameter estimate and its $1\sigma$ uncertainty.

The TESS photometry exhibits clear eclipses with differing depths as well as out-of-eclipse variability
on the orbital period. We fix the primary effective temperature according to the value derived in
Sect.~\ref{sec:primary}. Although the field covered by the TESS FFI surrounding the target contains
numerous Gaia sources, the brightest nearby target is more than five magnitudes fainter. Following this, we assumed that the light curve of V772~Cas\ is not diluted by any additional light source and fixed
the third light parameter to zero in our {\sc PHOEBE} model.

Preliminary binary model fits with variable eccentricity indicated that this parameter does not
exceed 0.02 and is statistically consistent with zero. Consequently, we fixed eccentricity in the
orbit to zero in the final analysis. This can be contrasted with a marginal eccentricity of $e=0.17\pm0.10$ inferred by \citet{gandet:2008} from low-precision photographic radial velocity measurements. Since the secondary eclipse was not detected by that study, the eccentricity could not have been constrained by the light curve.

Additionally, we apply a Gaussian prior on $v_1\sin i$, as taken
from Sect.~\ref{sec:primary}. For each sampled parameter combination, we calculate the corresponding
component surface gravities $\log g_1$ and $\log g_2$, and use these (along with the fixed $T_{\rm eff,1}$
and sampled $T_{\rm eff,2}$) to interpolate both limb-darkening and gravity-darkening coefficients for the
TESS pass-band from the tables provided by \citet{claret:2017}. For limb-darkening, we apply the square-root law.
Interpolating the gravity darkening coefficients according to the effective temperatures and surface
gravities of the components of a given model allows us to reduce the parameter space as opposed to varying
these coefficients and checking that they are theoretically consistent \textit{a posteriori}.

In addition to these fixed parameters, we vary the orbital period, $P_{\rm orb}$, the reference time, $T_0$, the
orbital inclination, $i$, the semi-major axis, $a$, the mass-ratio, $q$, the systemic velocity, $\gamma$,
as system parameters. We also vary potentials ($\Omega_{1,2}$), albedos ($A_{1,2}$), and light contributions ($l_{1,2}$) of the primary and secondary, and vary the secondary effective temperature as well as the primary synchronicity parameter. All these parameters are allowed to vary with uniform priors, with bounds (where applicable) listed in Table~\ref{tab:binary_pars}.

\begin{table}
    \centering

    \renewcommand{\arraystretch}{1.3}
        \caption{Sampled and derived binary parameters with boundaries and estimates
                 according to median values and uncertainties listed as 68 per cent HPD intervals.}

    \begin{tabular}{l c c l}
    \hline\hline
    Parameter &  & Prior range & HPD estimate \\
    \hline
    $P_{\rm orb}$ & d & $\mathcal{U}$(0,-) & $5.0138 \pm 0.0001 $ \\
    $T_0$ & d & $\mathcal{U}$(-,-) & $2458803.4016 \pm 0.0004$ \\
    $q$ & $\frac{M_2}{M_1}$ & $\mathcal{U}$(0.01,1) &  $0.2343 \pm 0.0006$ \\
    $a$ & $R_{\odot}$ & $\mathcal{U}$(1,50) & $20.60 \pm 0.03 $ \\
    $\gamma$ & ${\rm km\,s^{-1}}$ & $\mathcal{U}$($-50$,50) & $-1.6 \pm 0.1  $ \\
    $i$ & $\deg$ & $\mathcal{U}$(45,90) & $ 85.40 \pm 0.3$ \\
    $T_{\rm eff,1}$ & ${\rm K}$ & N/A & $13\,800$ \\
    $T_{\rm eff,2}$ & ${\rm K}$ & $\mathcal{U}$(3500,50000) & $5750 \pm 150$ \\
    $A_1$ & & $\mathcal{U}$(0,1) & $0.5\pm0.3$ \\
    $A_2$ & & $\mathcal{U}$(0,1) & $0.4\pm0.03$ \\
    $\Omega_1$ &  & $\mathcal{U}$(3,20) & $4.50 \pm 0.03$\\
    $\Omega_2$ &  & $\mathcal{U}$(3,20) & $6.90 \pm 0.04$\\
    $\omega_{\rm rot,1}/\omega_{\rm orb}$ & & $\mathcal{U}$(0.1,3) & $0.459 \pm 0.007$ \\
    $l_1$ & per cent & $\mathcal{U}$(50,100) & $99.57 \pm 0.01$\\
    $l_2$ & per cent & $\mathcal{U}$(0,50) & $0.43 \pm 0.01$\\
    \hline
    $r_1/a$ & & N/A & $0.234 \pm 0.001$ \\
    $M_1$ & $M_{\odot}$ & N/A & $3.78\pm0.02$ \\
    $R_1$ & $R_{\odot}$ & N/A & $4.83\pm0.03$ \\
    $\log g_1$ & cm\,s$^{-2}$ & N/A & $3.648\pm0.005$ \\
    $r_2/a$ & & N/A & $0.0424 \pm 0.0003$ \\
    $M_2$ & $M_{\odot}$ & N/A & $0.89\pm0.01$ \\
    $R_2$ & $R_{\odot}$ & N/A & $0.87\pm0.01$ \\
    $\log g_2$ & cm\,s$^{-2}$ & N/A & $4.50\pm0.01$ \\

    \hline
    \end{tabular}
    \label{tab:binary_pars}
\end{table}

\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{v772Cas_lc_rv_model.pdf}
\caption{{\bf Top}: TESS FFI observations (black symbols) and optimized {\sc PHOEBE} light curve model (red line) as a function of the orbital phase. {\bf Bottom}: Radial velocities derived from HERMES spectra (black symbols) and optimized {\sc PHOEBE} radial velocity model (red line).}
\label{fig:eb_soln}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{v772Cas_ecls}
\caption{{\bf Top}: TESS FFI observations (black symbols) and optimized {\sc PHOEBE} light curve model (red line) in the vicinity of primary (left) and secondary (right) eclipses. {\bf Bottom}: Residuals from the light curve fit.}
\label{fig:eclp}
\end{figure}

\begin{figure*}
\centering
\figps{\hsize}{90}{SED}
\caption{Comparison between the observed (symbols) and theoretical (lines) spectral energy distribution of V772~Cas\ in the ultraviolet and optical (left panel) and near-infrared (right panel) wavelength regions. Calculation shown with the thick solid line corresponds to primary's parameters $T_{\rm eff}$\,=\,13800~K, $\log{g}$\,=\,3.7, $R=4.95$~$R_\odot$, and reddening $E(B-V)=0.144$. Dotted lines illustrate the impact of changing $T_{\rm eff}$\ by $\pm500$~K.}
\label{fig:sed}
\end{figure*}

Despite the lack of detection of the secondary in the spectra or LSD profiles, the presence of
flat bottom secondary eclipse (which indicates a total eclipse) enables us to obtain some extra constraints on the stellar radii, and through the
determination of the potentials, the mass ratio. Additionally, the presence of ellipsoidal
variability provides further constraint on the mass ratio. The optimised parameter estimates and
their highest posterior density (HPD) 68 per cent ($1\sigma$) uncertainties are listed in Table~\ref{tab:binary_pars}.
The posterior distributions are illustrated in Figs.~\ref{fig:pd1} to \ref{fig:pd4}.
Our solution reports an evolving
intermediate mass $M_1=3.78\pm0.02$~$M_{\odot}$, $R_1=4.83\pm0.03$~$R_{\odot}$ primary and
a lower mass $M_2=0.89\pm0.01$~$M_{\odot}$, $R_2=0.87\pm0.01$~$R_{\odot}$ secondary which is
still in the first half of its main-sequence evolution. The size of both components is well below their Roche radii (10.8~$R_{\odot}$ for the primary and 5.0~$R_{\odot}$ for the secondary, respectively).

The optimized light curve and radial velocity models constructed from these parameters are shown in the top and bottom panels of Fig.~\ref{fig:eb_soln}, respectively. This figure demonstrates that the PHOEBE binary system model successfully reproduces available photometric and spectroscopic observations, both within the primary and secondary eclipses and outside the eclipses. The latter 0.6 per cent peak-to-peak photometric variation is thus interpreted as an ellipsoidal variability caused by a slightly distorted shape of the primary. Our analysis indicates that its maximum deviation from a spherical shape is about 0.7 per cent. In addition, Fig.~\ref{fig:eclp} shows observations, the light curve model, and the corresponding residuals around the primary and secondary eclipses. The residuals shown in the lower panel of Fig.~\ref{fig:eclp} reveal no systematic trends.

We do not include spots
in our binary model as adding those would introduce several additional free parameters. There is, in fact, no evidence of spot signatures in the photometric data after ellipsoidal variation is accounted for. Inspection of the residuals after subtraction of the binary light curve model does not reveal any clear indication of rotational modulation or any other periodic variability above $\approx100$~ppm for the frequencies below 2 d$^{-1}$ and 20~ppm in the 2--24 d$^{-1}$ frequency range.

Given the sharp points of ingress and egress (the first and fourth contacts) at both primary and secondary eclipse, we ran two {\sc PHOEBE} models with different fine and coarse grid sizes to characterise the numerical noise in our solution. We find a difference of 18~ppm between the two models, whereas the residual scatter is 230~ppm, meaning that the numerical noise is responsible for less than 10 per cent of the residual scatter \citep{maxted:2020}.

\begin{figure}
\centering
\figps{0.85\hsize}{0}{Hlines}
\caption{Comparison between the observed (symbols) and theoretical (lines) hydrogen Balmer line profiles. Calculations shown with solid lines employ $T_{\rm eff}$\,=\,13800~K, $\log{g}$\,=\,3.7 model atmosphere for the primary together with the chemical abundances determined in this study. Dotted lines illustrate the effect of changing $\log{g}$\ by $\pm0.2$~dex.}
\label{fig:hlines}
\end{figure}

\begin{figure*}
\centering
\figps{\hsize}{0}{spectral_fits}
\caption{Comparison of the average spectrum of V772~Cas\ (symbols) with the best fitting theoretical model spectrum (solid line) in several wavelength regions containing spectral lines commonly enhanced in HgMn stars. Dotted line shows synthetic spectrum calculated with solar abundances.}
\label{fig:fits}
\end{figure*}

\begin{figure}
\centering
\figps{\hsize}{0}{abundances}
\caption{Abundances of individual elements in the primary component of V772~Cas\ (filled symbols) relative to the solar chemical composition. Filled circles correspond to neutral and singly ionised species; the filled square shows the abundance of \ion{Si}{iii}. Open symbols illustrate abundances of HgMn stars with $T_{\rm eff}$\,=\,13000--14000~K from the catalogue by \citet{ghazaryan:2018}.}
\label{fig:abund}
\end{figure}

\section{Analysis of HgMn primary}
\label{sec:primary}

\subsection{Atmospheric parameters}

Here we use a combination of modelling the stellar spectral energy distribution (SED) and the hydrogen Balmer lines to determine primary's $T_{\rm eff}$\ and $\log{g}$\ respectively. This approach is commonly employed for normal and peculiar late-B and A-type stars (e.g. \citealt*{ryabchikova:1999a}; \citealt{fossati:2009}; \citealt{rusomarov:2016}; \citealt*{kochukhov:2019}). Other spectral indicators, such as \ion{He}{i} and metal lines, cannot be used for the atmospheric parameter determination of HgMn stars due to non-solar photospheric element abundances and occasional vertical chemical stratification.

Model atmosphere analysis of the primary component of V772~Cas\ was carried out using the {\sc LLmodels} code \citep{shulyak:2004}, taking into account individual atmospheric abundances. The influence of the faint secondary was ignored. The effective temperature was determined by comparing the model SED with the TD1 satellite stellar flux measurements in UV \citep{thompson:1978} as well as the optical and near-infrared fluxes obtained from Geneva \citep{hauck:1982} and 2MASS \citep{cutri:2003} photometric measurements, respectively. The reddening $E(B-V)=0.144\pm0.050$ \citep{lallement:2019} and the Gaia DR2 distance $363\pm6$~pc \citep{gaia-collaboration:2018} were adopted for SED fitting with the stellar effective temperature and radius adjusted to match the observations. This analysis yielded $T_{\rm eff}$\,=\,$13800\pm500$~K and $R=4.95\pm0.2$~$R_\odot$. Figure~\ref{fig:sed} compares observations with the theoretical SED computed with {\sc LLmodels} using these parameters. The stellar effective temperature derived from the SED agrees reasonably well with $T_{\rm eff}$\,=\,13400--13800~K that can be obtained for this star using different Str\"omgren and Geneva photometric calibrations (\citealt{kunzli:1997}; \citealt*{paunzen:2005a}; \citealt*{paunzen:2006}).

The surface gravity was determined by fitting the observed hydrogen Balmer line profiles in the time-averaged HERMES spectrum with the theoretical calculations using the {\sc Synth3} spectrum synthesis code \citep{kochukhov:2007d} and the {\sc LLmodels} atmospheres described above. The average spectrum was constructed by co-adding ten individual observations after correcting the radial velocity shifts reported in Table~\ref{tab:obs}. As demonstrated by Fig.~\ref{fig:hlines}, the hydrogen lines in the mean spectrum of the primary are well approximated with $\log{g}$\,=\,$3.7\pm0.1$. This spectroscopic estimate of $\log{g}$\ and the stellar radius inferred from the SED are consistent within uncertainty with the results of binary system modelling with {\sc PHOEBE} in Sect.~\ref{sec:binary_results}.

\subsection{Abundances}

Chemical abundances were estimated by fitting {\sc Synth3} spectra to short segments of the average spectrum of V772~Cas. In these fits, individual abundances of one or several elements as well as the projected rotational velocity were determined using the {\sc BinMag} \citep{kochukhov:2018}\footnote{\url{https://www.astro.uu.se/~oleg/binmag.html}} IDL tool. Different lines of the same ions were analysed independently and the scatter of abundances, quantified by the standard deviation, was taken as an uncertainty estimate. The input line list for these calculations was extracted from VALD, using the latest version of the database that incorporates hyperfine and isotopic splitting \citep{pakhomov:2019}. For V772~Cas\ the former is particularly important for accurate analysis of \ion{Mn}{ii} and \ion{Ga}{ii} lines. Information on the hyperfine and isotopic splitting of the \ion{Hg}{ii} 3984~\AA\ line was taken from \citet{woolf:1999}.

Several examples of spectrum synthesis fits around the lines typically enhanced in HgMn stars are shown in Fig.~\ref{fig:fits}. Examples of fits to other ions analysed in the paper are demonstrated in Fig.~\ref{fig:other_fits}. We were able to obtain abundance estimates for 19 ions based on 164 individual lines and blends, listed in Table~\ref{tab:lines}. The resulting abundances are presented in Table~\ref{tab:abund}, which lists the ions studied, the number of lines analysed, the average abundance and corresponding error for V772~Cas, the corresponding average abundance of HgMn stars with $T_{\rm eff}$\,=\,13000--14000~K from the compilation by \citet{ghazaryan:2018}, and the solar abundance of each ion \citep{asplund:2009}. Our abundance analysis of V772~Cas\ was carried out under the local thermodynamic equilibrium (LTE) assumption. Several abundances reported in Table~\ref{tab:abund}, for example Ca \citep*{sitnova:2018}, Si \citep{mashonkina:2020}, and Ne \citep{alexeeva:2020}, are likely affected by departures from LTE leading to abundance corrections of about 0.1--0.3~dex.

The abundance pattern of V772~Cas\ relative to the solar chemical composition is illustrated in Fig.~\ref{fig:abund}. We also show in this figure abundances of 33 HgMn stars in the 13000--14000~K $T_{\rm eff}$\ range from the catalogue by \citet{ghazaryan:2018}. It is evident that V772~Cas\ exhibits an unremarkable abundance pattern, very often seen in HgMn stars with a similar $T_{\rm eff}$. In particular, He is underabundant, Si is close to solar, P and Mn are overabundant by $\sim$\,2~dex, and Ga, Xe, and Hg are overabundant by up to 4~dex. Pr is also overabundant in V772~Cas\ by $\approx$\,2.5~dex, as found for several other HgMn stars, although its abundance estimate based on a single blended \ion{Pr}{iii} line is somewhat uncertain.

Considering the projected rotational velocity, determined with the spectrum synthesis fit from 84 unblended lines, we found $v_{\rm e} \sin i$\,=\,$22.3\pm0.3$~km\,s$^{-1}$. Together with the stellar radius $R=4.83\pm0.03$~$R_\odot$\ and inclination of the rotation axis $i_{\rm rot}=i_{\rm orb}=85.4\pm0.3$\degr\ found in Sect.~\ref{sec:binary}, this implies $P_{\rm rot}=10.9\pm0.2$~d. Thus, the primary is rotating sub-synchronously with a period about twice longer than the orbital period provided that the orbital and rotational axes are aligned. The latter assumption is reasonable for close binaries on theoretical grounds \citep{zahn:1977,hut:1981} and generally agrees with observational findings \citep{hale:1994,farbiash:2004}.

\begin{table}
\caption{Atmospheric chemical composition of the primary component of V772~Cas. The columns give the ion identification, the number of lines studied, abundance for V772~Cas, average abundance for HgMn stars in the 13000--14000 $T_{\rm eff}$\ interval \citep{ghazaryan:2018}, and the solar abundance \citep{asplund:2009}. \label{tab:abund}}
 \begin{tabular}{lllll}
\hline\hline
Ion & $N$ & \multicolumn{3}{c}{$\log(N_{\rm el}/N_{\rm tot})$} \\
\cline{3-5} \\[-0.2cm]
& & V772~Cas\ & $\langle \mathrm{HgMn} \rangle$ & Sun \\
\hline
\ion{He}{i}   &   5 & $-2.02\pm0.10$  & $-2.01\pm0.52$ & $ -1.11$ \\
\ion{C}{ii}   &   3 & $-3.80\pm0.15$  & $-4.05\pm0.48$ & $ -3.61$ \\
\ion{O}{i}    &   2 & $-3.37\pm0.16$  & $-3.36\pm0.13$ & $ -3.35$ \\
\ion{Ne}{i}   &  11 & $-3.63\pm0.04$  & $-4.29\pm0.47$ & $ -4.11$ \\
\ion{Mg}{ii}  &   3 & $-4.73\pm0.14$  & $-5.16\pm0.53$ & $ -4.44$ \\
\ion{Si}{ii}  &   5 & $-4.38\pm0.09$  & $-4.89\pm0.63$ & $ -4.53$ \\
\ion{Si}{iii} &   2 & $-4.56\pm0.01$  & $-4.89\pm0.63$ & $ -4.53$ \\
\ion{P}{ii}   &  21 & $-4.63\pm0.22$  & $-4.94\pm0.48$ & $ -6.63$ \\
\ion{S}{ii}   &   8 & $-5.89\pm0.10$  & $-5.28\pm0.63$ & $ -4.92$ \\
\ion{Ca}{ii}  &   1 & $-5.10       $  & $-5.32\pm0.41$ & $ -5.70$ \\
\ion{Ti}{ii}  &   8 & $-6.10\pm0.07$  & $-6.58\pm0.36$ & $ -7.09$ \\
\ion{Cr}{ii}  &   6 & $-5.68\pm0.06$  & $-6.24\pm0.51$ & $ -6.40$ \\
\ion{Mn}{ii}  &  38 & $-4.17\pm0.14$  & $-4.97\pm0.86$ & $ -6.61$ \\
\ion{Fe}{ii}  &  42 & $-4.11\pm0.10$  & $-4.32\pm0.50$ & $ -4.54$ \\
\ion{Ni}{ii}  &   1 & $-5.91       $  & $-6.20\pm0.43$ & $ -5.82$ \\
\ion{Ga}{ii}  &   3 & $-5.42\pm0.20$  & $-5.87\pm1.13$ & $ -9.00$ \\
\ion{Xe}{ii}  &   3 & $-5.37\pm0.35$  & $-5.51\pm0.58$ & $ -9.80$ \\
\ion{Pr}{iii} &   1 & $-8.66       $: & $-9.21\pm1.10$ & $-11.32$ \\
\ion{Hg}{ii}  &   1 & $-6.81       $  & $-6.49\pm1.09$ & $-10.87$ \\
\hline
\end{tabular}
\end{table}

\section{Discussion}
\label{sec:disc}

In this paper we investigated the nature of the bright but poorly studied eclipsing binary system V772~Cas. Based on the new high-precision photometric data provided by the TESS mission we confirmed the presence of primary eclipses and identified secondary eclipses for the first time. Significant out-of-eclipse photometric variability, synchronised with the orbital motion, was confirmed for this system. We have acquired high-resolution spectra of V772~Cas\ with the aim to determine fundamental parameters of the components and better characterise surface abundance pattern of the primary star, which in the past was attributed conflicting HgMn and BpSi spectral classifications.

Our analysis reveals that V772~Cas\ is an SB1 system with a late-B primary showing abundance anomalies typical of HgMn stars. Measurement of the projected rotational velocity shows that the primary rotates sub-synchronously. Taking into account constraints from the atmospheric modelling of the primary, we carried out detailed modelling of the TESS light curve with the PHOEBE code. This analysis demonstrated that the conspicuous out-of-eclipse variability is explained by the ellipsoidal shape of the primary, which arises due to the tidal interaction with the secondary. The ellipsoidal variability is clearly visible in the TESS light curve despite a relatively low mass of the secondary thanks to a high precision of the space photometry data. Interpretation of the ellipsoidal variability together with eclipses and radial velocity variation of the primary allowed us to retrieve accurate masses and radii of both components despite the absence of the secondary's lines in the optical spectra. We concluded that an evolved HgMn primary of V772~Cas\ is orbited by a main sequence solar-type secondary, resulting in large luminosity and mass ratios.

Our results confirm that V772~Cas\ should be formally classified as an Algol-type (EA) eclipsing binary \citep{samus:2017}, without any additional types of variability present. In particular, the absence of rotational modulation due to surface spots and the HgMn spectroscopic classification of the primary demonstrates that this star is not an $\alpha^2$~CVn-type magnetic variable similar to the magnetic Bp stars recently identified in the eclipsing binaries HD~66051 \citep{kochukhov:2018b} and HD~65658 \citep{shultz:2019}. Instead, V772~Cas\ should be discussed in the context of research on HgMn stars. These objects lack strong large-scale magnetic fields \citep{shorlin:2002,auriere:2010a,bagnulo:2012,kochukhov:2011b,makaganiuk:2011a,makaganiuk:2011,makaganiuk:2012} and also possess no complex tangled fields stronger than a few hundred G \citep{kochukhov:2013a}. Despite this, these stars are able to develop a low-contrast non-uniform surface abundance distribution of some heavy elements \citep{adelman:2002,kochukhov:2005,folsom:2010,makaganiuk:2011}. Geometry of these abundance spots appears to slowly evolve with time \citep{kochukhov:2007b,briquet:2010,korhonen:2013}. Photometric time series studies of HgMn stars occasionally reveal rotational modulation, presumably related to heavy-element abundance spots (\citealt{morel:2014}; \citealt{strassmeier:2017}; \citealt{white:2017}; \citealt*{prvak:2020}), and, possibly, to SPB pulsations \citep{hummerich:2018}. Our study shows that ellipsoidal variability is yet another phenomenon that can contribute to or even dominate the light curves of HgMn stars in close binaries.

\begin{table*}
\caption{Properties of eclipsing binary systems containing HgMn stars. When available, uncertainties in the last significant digit are indicated by the numbers in brackets. \label{tab:hgmns}}
 \begin{tabular}{lllllllllll}
\hline
\hline
Star & $V$ & Spec & $P_{\rm orb}$ (d) & $e$ & $q$ & $M_1$ ($M_\odot$) & $R_1$ ($R_\odot$) & $M_2$ ($M_\odot$) & $R_2$ ($R_\odot$) & Reference \\
\hline
AR Aur          &   6.14 & SB2 &    4.13 & 0.0    & 0.927(2) & 2.544(9) & 1.80(1) & 2.358(8) & 1.83(2) &  \citet{hubrig:2012} \\
TYC 455-791-1   &  11.95 & SB2 &   12.47 & 0.18   & 0.941(8)  & 3.1  & 2.4  & 2.9  & 2.3  &  \citet{strassmeier:2017} \\
V772 Cas        &   6.68 & SB1 &    5.01 & 0.0    & 0.2343(6) & 3.78(2) & 4.83(3) & 0.89(1) & 0.87(1) &  This work \\
\hline
\end{tabular}
\end{table*}

Although HgMn stars are commonly found in close binary systems \citep{abt:1973,gerbaldi:1985}, only two eclipsing binaries containing HgMn components were known prior to this study. They are AR~Aur \citep{nordstrom:1994,folsom:2010,hubrig:2012} and TYC 455-791-1 \citep{strassmeier:2017}. Table~\ref{tab:hgmns} summarises the properties of these three systems in comparison to V772~Cas. All three eclipsing binaries are non-interacting, detached systems, so their components should have evolved as if they were single stars (\citealt*{torres:2010}). The resulting mass-radius relationship is compared to single-star MESA (Modules for Experiments in Stellar Astrophysics) theoretical stellar evolution models from the MIST grid\footnote{\url{http://waps.cfa.harvard.edu/MIST}} \citep{dotter:2016,choi:2016} in Fig.~\ref{fig:isoc}. These models are available for two values of initial rotational velocity, a wide metallicity range, and a single empirically calibrated overshooting prescription, consistent with recent observational constraints \citep{claret:2019}. In this work, we have chosen to use the isochrones that were computed ignoring the stellar rotation. This assumption is appropriate for the three systems considered here and for HgMn stars in general since a slow rotation is known to be a necessary condition for this chemical peculiarity to appear \citep{michaud:1982}.

The chemical anomalies produced by the atomic diffusion in A and B stars are constrained to the outermost stellar layers and are not indicative of the bulk metallicity (e.g. \citealt*{richard:2001}). The latter is unknown but is believed to not differ much from that of normal stars. According to \citet{sofia:2001}, the metallicity scatter of young F and G disk stars in the solar neighbourhood reaches 0.1--0.2~dex, which is similar to the 0.4~dex metallicity range of open clusters containing ApBp stars \citep{bagnulo:2006,landstreet:2007}. Therefore, we included a metallicity variation by $\pm0.2$~dex around the solar value when estimating the age of V772~Cas\ from the MIST isochrones.

In all three systems the star with HgMn peculiarity is the primary. AR~Aur and TYC 455-791-1 are close to ZAMS and have mass ratios not far from unity. In fact, AR~Aur~B is probably still contracting towards ZAMS, indicating the extreme youth of this system \citep{nordstrom:1994}. The masses of the components of AR~Aur and TYC 455-791-1 span a narrow range of 2.4--3.1 $M_\odot$. On the other hand, the primary of V772~Cas\ is more massive and considerably more evolved. It is likely to be at or near the TAMS whereas the solar-type secondary is still on the main sequence. The location of V772~Cas\ components on the mass-radius diagram illustrated in Fig.~\ref{fig:isoc} is best described by a set of isochrones with ages from $\log t/\mathrm{yr}=8.26$ to 8.32 for the metallicity range discussed above and the standard overshooting prescription adopted in the MIST grid. A similar age range (8.32--8.38) is obtained by applying the isochrone-cloud fitting methodology \citep{johnston:2019} with a different grid of solar-metallicity MESA models spanning a large range of overshooting parameter.

\begin{figure}
\centering
\figps{\hsize}{0}{hgmn_isoc}
\caption{Masses and radii of the eclipsing binaries with HgMn components. Parameters of the primaries are indicated by filled symbols. Open symbols correspond to the secondaries. The three systems included in this plot are V772~Cas\ (circles), AR~Aur (squares), and TYC 455-791-1 (diamonds). Theoretical MESA isochrones corresponding to ages from $\log t/\mathrm{yr}=7.5$ to 8.8 are also shown. The underlying grey curves illustrate the effect of varying metallicity by $\pm0.2$~dex.}
\label{fig:isoc}
\end{figure}

Despite this dramatic difference in evolutionary status, the HgMn primaries of all three eclipsing binary systems share the same basic chemical abundance pattern. This observation suggests that the HgMn-type chemical peculiarity develops rapidly at or near the ZAMS and then persists throughout the entire main sequence life time of an intermediate-mass star. This is broadly in agreement with the investigation of H-R diagram positions of single HgMn stars by \citet*{adelman:2003a}, although that study showed that these stars tend to concentrate in the first half of the main-sequence life time and that there are very few luminous HgMn stars comparable to V772~Cas\ A.

\section*{Acknowledgements}
Based on observations made with the Mercator Telescope, operated on the Island of La Palma by the Flemish Community, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrof\'isica de Canarias.
This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST).
O.K. acknowledges support by the Swedish Research Council and the Swedish National Space Board.
J.L.-B. acknowledges support from FAPESP (grant 2017/23731-1).
This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by NASA's Science Mission directorate.
This research made use of Astropy,\footnote{\url{http://www.astropy.org}} a community-developed core Python package for Astronomy \citep{astropy:2013,astropy:2018}
as well as the Corner\footnote{\url{https://corner.readthedocs.io}} Python code by \citet{foreman-mackey:2016}.
The research leading to these results has received funding from the European
Research Council (ERC) under the European Union's Horizon 2020 research and
innovation programme (grant agreement N$^\circ$670519: MAMSIE), from the
KU\,Leuven Research Council (grant C16/18/005: PARADISE), from the Research
Foundation Flanders (FWO) under grant agreements G0H5416N (ERC Runner Up Project)
and G0A2917N (BlackGEM), as well as from the BELgian federal Science Policy Office
(BELSPO) through PRODEX grant PLATO.
T.R. thanks the Ministry of Science and Higher Education of Russian Federation
(grant 13.1902.21.0039) for partial financial support.
M.E.S. acknowledges  support from the Annie Jump Cannon Fellowship, supported by the University of Delaware and endowed by the Mount Cuba Astronomical Observatory.

\section*{Data availability}
The data underlying this article will be shared on reasonable request to the corresponding author.


\section{Introduction}

The spread of COVID-19 pandemic has practically affected the entire planet, and created enormous challenges on every aspect of human life and organization, starting with the health sector and with far reaching consequences on the economy, education, sports, transportation and politics. Since the first cases were registered in Wuhan, China in December 2019 \cite{WHO}, the global spatial dynamics of the infection have been changing as the disease swiftly moved towards the West \cite{Mathys1} into Europe then into the United States, South America, and eventually to the whole world, with nearly 38.1 million cases and 1.1 million deaths registered until October 12, 2020 \cite{Worldometer}.

Given the global geographic spread of the virus and the local wide spread in many countries, and the nature of the transmission of the virus, it is important to understand the spatial mechanisms of this spread and its dependence on proximity, demographics and social characteristics of infected areas. Spatial analysis provides a better understanding of the routes of transmission of infections \cite{Meng1}, consequently, it allows the decision-makers to draft and implement effective health and mitigation measures to reduce risks associated with the pandemic.

In Lebanon, the first case was registered on February 21, 2020 \cite{MOPH} and by October 12, $54624$ cases and $466$ deaths were registered \cite{DRM}. The first few weeks witnessed a relatively rapid increase but it sharply declined as a result of the strong mitigation measures enforced by the beginning of March. The lift of the international travel ban and the partial easing of measures led to the revival of higher spread rates since July. Only $1788$ cases were registered by July 1, 2020 before a sharp rise from July through October. The cases were mainly concentrated in Beirut, its suburbs and its neighboring areas in Mount Lebanon. On August 4, a huge explosion rattled the port of Beirut and destroyed thousands of houses and buildings in the surrounding areas. People were rushed into hospitals, with thousands of injuries recorded on that day \cite{BBC}. On such a horrible incident, hundreds of volunteers and civil defense teams were involved in rescue work for several days. The social distancing measures were largely neglected in such an emergency situation. The spread accelerated in the upcoming weeks, with sharp rise in Beirut and its surroundings and with a national widespread reaching all regions and major towns and cities \cite{Science}.

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.8]{Infection-spread}
\caption{A monthly map of the regional cumulative number of infections in Lebanese districts between April and October 2020.}\label{fig:1}
\end{center}
\end{figure}

\emph{Related Literature}: Spatial autocorrelation is the statistical analysis of data studied in space or in space-time aiming for the identification and estimation of spatial processes \cite{CliffOrd1, Anselin}. It has been implemented to study and analyze the spread of various diseases and infections including cancer, diabetes, SARS, influenza virus, COVID-19, etc... \cite{Glick, Hipp, SARS, PLOS, Kuwait}. Recent studies also inspected the effect of city size, population, transportation systems and demographics on the disease spread and its mortality rate \cite{PLOS1, Urban, Xie, Safety1, Safety2}.   The determination of the mean center of a population (centroid) was discussed in \cite{Aboufadel, Bachi, Barmore} and extending the concept to the determination of the mean center of wealth and infections allowed for a spatial analysis of the temporal dynamics of wealth distribution, economic growth and infectious diseases \cite{Mathys2}. The dynamics of the outbreak of COVID-19 in Lebanon and its reproduction number dynamics were studied in \cite{DeebJalloul1, DeebJalloul2, Mourad}.

In this paper, we study the clustering and spatial progression of new infections in Lebanon by applying the methods of spatial autocorrelation with different model parameterizations of geographic, demographic and social variables including  adjacency, proximity, population, population density, poverty rate and poverty density. Locating the mean center of the epidemic spread as a function of time is used to analyze the temporal geographic development of the spread. The obtained results provide a solid basis for the concerned policy makers to draw well-grounded and scientifically based local and regional measures that would contribute to controlling the infection spread.

 The paper is organized as follows: in section 2 we introduce the implemented analytic mathematical and statistical methods and tools. Results are presented and discussed in section 3, and section 4 concludes the paper.

\section{Analytic methods and tools}

\subsection{Moran's I index}

Moran's $I$ index is an inferential statistic used to measure the spatial autocorrelation based both on  locations and feature values simultaneously. It is defined as \cite{CliffOrd1}:
\begin{equation}
 I=\frac{N\Sigma_{ij}W_{ij}(X_{i}-\bar{X})(X_{j}-\bar{X})}{\Sigma_{ij}W_{ij}\Sigma_{i}(X_{i}-\bar{X})^{2}}
\end{equation}

where $W_{ij}$ represents different types of adjacency between
region $i$ and region $j$, corresponding to different models of
infectious spread. $N$ is the number of regions under consideration
and $X_{i}$ represents the number of new daily infections
in district $i$. $\bar{X}$ is the average number
of new daily infections per region, and it is given by $\bar{X}=\frac{\Sigma_{i}X_{i}}{N}$.
The numerical outcome of $I$ falls between $-1$ and $1$ and it indicates whether a distribution is dispersed, random or clustered. A value of $I$ close to $0$ indicates a random distribution, while positive values  indicate clustered spatial distribution and negative values indicate dispersion. Larger values of $|I|$ nearer to $1$ mean stronger clustering (positive $I$) or stronger dispersion (negative $I$).

The $z_{I}$-score associated to this statistic is defined by:

\begin{equation}
z_{I}=\frac{I-E[I]}{\sqrt{V[I]}}
\end{equation}

where the expected value $E[I]$ and the variance $V[I]$ are defined in the Appendix. The $z$-score or the corresponding $p$-value of the statistic are used to reject the null hypothesis and eliminate the possibility of a random pattern leading to the obtained value of the Moran $I$ statistic.

In this paper, we take a $95\%$ confidence level corresponding to $|z_{I}|>1.96$ or equivalently to $p<0.05$ in order to confirm the outcome of clustering or dispersion of our spatial data indicated by $I$. In this case we say that the $p$-value is statistically significant, and based on the value of $I$ we can determine the pattern of the distribution.

We consider a model with six different cases of parameterization of the adjacency matrix $W_{ij}$ corresponding to geographic adjacency (case I), proximity (case II), population (case III), population density (case IV), poverty rate (case V) and poverty density (case VI). Table \ref{table:1} summarizes relevant data from the Lebanese districts.

In  casel I, we take $W_{ij}=1$ for districts sharing common borders and $W_{ij}=0$ otherwise, while in case II we determine $W_{ij}=\frac{1}{d_{ij}}$ where $d_{ij}$ is the driving distance between the administrative centers of regions $i$ and $j$. Those two cases study the effect of administrative adjacency and the distance proximity of different districts on the geographic clustering of new infections in Lebanon.

In case III and case IV, we apply the methods used in (\cite{Meng1, Kang}) to analyze the effects of population and population density on the spread of the disease since the virus is carried by people and its spread is supposed to be related to their population and interaction. We sort the districts by the number of their residents (obtained from \cite{CAS})  and then by the density of their residents relative to their areas. In these two cases, districts of consecutive populations and population densities  are assigned a factor of $W_{ij}=1$, and $W_{ij}=0$ otherwise. This provides a statistic about the clustering of infections according to population and population density respectively.

Lastly, in cases V and VI, we introduce  new parameters, namely the poverty rate and the poverty density in different districts and we analyze their effect on infection clustering. We sort the distritcs by their rates of poverty  and poverty density  \cite{CAS} and assign  $W_{ij}=1$ for regions of consecutive order of poverty rate or poverty density, and $W_{ij}=0$ otherwise, in a similar methodology to cases III and IV in order to infer the effect of poverty rate and density on the geographical patterns of infection spread.

\subsection{Mean center of infection}

The mean center of infection (henceforth MCI) is a geographic location that represents the weighted mean of the positions of infected individuals on the surface of Earth, assumed to be spherical. Assigning the value of Earth's radius to unity, the two spherical coordinates that determine the unique position of a point are its latitude $\lambda_{i}$ and longitude $\phi_{i}$.  The latitude is a measurement of location north or south of the equator while the longitude is a measurement of location east or west of the prime meridian at Greenwich, UK.

The Cartesian position vector $\overrightarrow{r_{i}}=(x_{i},y_{i},z_{i})$
is related to spherical coordinates with unit radius by the relations \cite{Coordinates}:

\begin{equation}
\begin{cases}
x_{i}=\cos\lambda_{i}\cos\phi_{i}\\
y_{i}=\sin\lambda_{i}\cos\phi_{i}\\
z_{i}=\sin\phi_{i}
\end{cases}
\end{equation}

We denote the district number of infections (new or cumulative) by $X_{i}$ as defined above, and the Cartesian positions of the administrative centers by $(x_{i},y_{i},z_{i})$ . Then, the Cartesian position of the weighted mean of infections
$\overrightarrow{\hat{r_{i}}}$ is given by:

\begin{equation}
\begin{cases}
\hat{x}=\frac{\Sigma_{i}X_{i}x_{i}}{\Sigma_{i}X_{i}}\\
\hat{y}=\frac{\Sigma_{i}X_{i}y_{i}}{\Sigma_{i}X_{i}}\\
\hat{z}=\frac{\Sigma_{i}X_{i}z_{i}}{\Sigma_{i}X_{i}}
\end{cases}
\end{equation}

As suggested by \cite{Aboufadel}, the precise position on the surface of a sphere can
be determined from the normalized position vector defined by $\overrightarrow{\bar{r_{i}}}=(\bar{x},\bar{y},\bar{z})=\frac{\overrightarrow{\hat{r_{i}}}}{\mid\overrightarrow{\hat{r_{i}}}\mid},$leading
to

\begin{equation}
\begin{cases}
\bar{x}=\frac{\hat{x}}{\sqrt{\hat{x}^{2}+\hat{y}^{2}+\hat{z}^{2}}}\\
\bar{y}=\frac{\hat{y}}{\sqrt{\hat{x}^{2}+\hat{y}^{2}+\hat{z}^{2}}}\\
\bar{z}=\frac{\hat{z}}{\sqrt{\hat{x}^{2}+\hat{y}^{2}+\hat{z}^{2}}}
\end{cases}
\end{equation}

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.55]{Neighbourhood}
\caption{The figure shows Moran's I index and its corresponding p-value for cases I and II accounting for adjacency and proximity of new infections registered in Lebanese districts. The test shows that since July 2020, there is clear geographic clustering of infections. }\label{fig:2}
\end{center}
\end{figure}

\begin{figure}[h]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.55]{Demographics}
\caption{Moran's I index and its corresponding p-value for cases III and IV accounting for population and population density of different districts. Population density contributes for strong clustering since August 2020, while the districts' population effect is weaker. } \label{fig:3}
\end{center}
\end{figure}

Consequently, we can recover the spherical position of the mean center of infections by calculating the mean latitude and longitude as:

\begin{equation}
\begin{cases}
\bar{\phi}=\sin^{-1}\bar{z}\\
\bar{\lambda}=\tan^{-1}\frac{\bar{y}}{\bar{x}}
\end{cases}
\end{equation}

The latitude and the longitude can be located and plotted on maps and geographic information systems. We employ the spherical coordinates of geographic locations of the capitals of the $26$ administrative districts in Lebanon and the number of daily and cumulative infections in each region in order to determine the daily MCI accordingly. This provides a tool to analyze the temporal dynamics of the mean geographic spread of the disease.

\section{Results and discussions}

The determination of the Moran's $I$ index and its corresponding p-value for the effect of adjacency and proximity of cases I and II on the clustering of daily new infections of COVID-19 in Lebanon shown in Figure (\ref{fig:2}), leads to the conclusion that since July 2020, there is strong clustering of infections in regions sharing common borders and among nearer regions. There were only few days when new infections were not clustered in adjacent regions, and only one day where distance was not shown to be a detrimental effect in the spatial spread of new cases. The maximum value of Moran's $I$ reached $0.660$ for case I and $0.380$ for case II indicating a high level of geographic clustering of the disease spread since July. The infections before July had a high $p$-value, indicating a high probability for random geographic spread.

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.55]{Poverty}
\caption{Moran's I index and its corresponding p-value for cases V and VI accounting for regional poverty rate and poverty density. The poverty rate is not a decisive factor for spatial spread but poverty density contributes to spatial clustering since the end of August 2020. }\label{fig:4}
\end{center}
\end{figure}

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.55]{Centroid}
\caption{This figure is a plot of the latitude and longitude of the weighted geographic center of COVID-19  infections in Lebanon. The upper graphs represent the temporal progression of the cumulative number of infections while the lower graphs represent that of the new daily cases.}\label{fig:5}
\end{center}
\end{figure}

The results of the spatial spread dynamics in relation to population and population density adjacency as shown in Moran's $I$ and $p$-value of cases III and IV depicted in  Figure (\ref{fig:3}) reveal that the spread was not clustered with respect to the regional population until late August 2020, where it started achieving a positive value of $I$ with $p<0.05$ indicating spatial clustering between regions of adjacent population rank, with several days showing a probability of random spread. The maximum attained $I$ was $0.666$. However, the statistics for districts with adjacent rank of population density show very strong spatial clustering since the middle of August with $I$ attaining a maximum value of $0.832$, which is the highest among all six studied cases.

The results of case V  (Figure \ref{fig:4}) show that the spatial spread cannot be attributed to adjacent ranking of poverty rates among the districts since the p-values remain above the $5\%$ level of confidence up until October 2020, hence no spatial clustering occurs. But when we consider the poverty density in case VI, we obtain positive values for Moran's $I$ since the end of August, with $p<0.05$ except for five days. Hence, spatial clustering among regions with adjacent ranking of poverty density occurs. The maximum attained $I$ in this case is $0.666$.

In comparison, we find out that clustering of new infections occurs starting on different dates between July and August for all considered cases except for case V corresponding to district populations. The strongest level of spatial clustering (highest $I$) occurs for model IV of population density after mid-August, while clustering associated to geographic adjacency and proximity (cases I and II) has the longest time span (since early July) and the highest levels of confidence.

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.65]{MCI}
\caption{The figure shows the variation of the geographic location of the MCI on a monthly basis between April 1 and October 1, 2020, together with the mean center of the population.}\label{fig:6}
\end{center}
\end{figure}

The location of the MCI was determined as a function of time as shown in   Figure (\ref{fig:5}). The mean latitude and longitude of the infection were determined according to the methods described in equation (6). The location of the cumulative MCI is plotted on the geographic map of Lebanon during the same period in  Figure (\ref{fig:6}), together with the mean center of population of the country. It started near the city of Jounieh, north-west of the mean center of population, but it has moved southward since May through August, where it started moving northward again. The location of the MCI of new infections was quite geographically distributed before July as the lower plot of  Figure (\ref{fig:5}) shows, before becoming more homogenous afterwards.

The reproduction number $R$ and the rate of the infection spread correlate with people's mobility \cite{Safety3}. Geographic clustering occurs because people's motion and local travel is higher in their close neighborhoods, especially in a country like Lebanon where with the absence of national public transportation throughout the country \cite{Transports} diminishes nationwide mobility. Higher levels of social interaction among people in dense regions also contribute to the spread of the disease, and this has shown the strongest clustering effect. The study of the spread of the infection allows the relevant authorities to draw appropriate country-specific and regional measures to curb the spread.

\begin{flushleft}

\begin{table}[t]
\begin{raggedright}
{\small{}
\caption{The table shows the distribution of the cumulative number of cases among the 26 Lebanese districts on October 12, 2020, with their respective populations, population densities, poverty rates and poverty densities.} \label{table:1}
\smallskip
\begin{tabular}{lc lc lc lc lc lc}
\hline
{\small{}Region Name} & {\small{}Number of cases} & {\small{}Population } & {\small{}Population denisty} & {\small{}Poverty} & {\small{}Poverty density}\tabularnewline
{\small{}} & {\small{} (Cumulative)} & {\small{} $\times100$} & {\small{}(resident/$km^{2}$)} & {\small{}rate (\%)} & {\small{}(resident/$km^{2}$)}\tabularnewline
\hline
{\small{}Akkar} & 1171 & 3204 & 418  & 38.4 & 38 \tabularnewline
{\small{}Minieh-Denniyeh} & 723  & 1408  & 389  & 48.6  & 189 \tabularnewline
{\small{}Tripoli} & 4198 & 2438 & 9030 & 31.7 & 2862 \tabularnewline
{\small{}Zgharta} & 854 & 877 & 399  & 25  & 100 \tabularnewline
{\small{}Koura} & 556 & 846 & 489 & 14.3 & 70 \tabularnewline
{\small{}Bcharre} & 115 & 221  & 140  & 13.4  & 19 \tabularnewline
{\small{}Batroun} & 341 & 589 & 212 & 5.5 & 12 \tabularnewline
{\small{}Jbeil} & 954 & 1295 & 301  & 12.7  & 38 \tabularnewline
{\small{}Kesrwan} & 1978 & 2605  & 762 & 18.5 & 141 \tabularnewline
{\small{}Meten} & 6139  & 5110  & 1928  & 20.4  & 393 \tabularnewline
{\small{}Beirut} & 6443  & 3417  & 17258  & 25.5 & 4401 \tabularnewline
{\small{}Baabda} & 7277  & 5538  & 2855 & 26.8  & 765 \tabularnewline
{\small{}Aley} & 3047 & 3008 & 1144  & 29.4  & 336 \tabularnewline
{\small{}Chouf} & 1965  & 2770  & 560  & 24.2  & 135 \tabularnewline
{\small{}Jezzine} & 125 & 321  & 133  & 21.9  & 29 \tabularnewline
{\small{}Saida} & 2472  & 2966  & 1079  & 19.1  & 206 \tabularnewline
{\small{}Nabatieh} & 683  & 1802  & 593  & 28.2  & 167 \tabularnewline
{\small{}Sour} & 1023 &  2557 & 933  & 30.3 & 283 \tabularnewline
{\small{}Bent Jbeil} & 331 & 962  & 364  & 22.9  & 83 \tabularnewline
{\small{}Marjeyoun} & 214 & 740  & 279  & 24.2  & 68 \tabularnewline
{\small{}Hasbaya} & 76 & 287  & 108  & 23.9  & 26 \tabularnewline
{\small{}Rachaya} & 79 & 338  & 62  & 16  & 10 \tabularnewline
{\small{}West Beqaa} & 474  & 864 & 184  & 25.5  & 47 \tabularnewline
{\small{}Zahleh} & 2367 & 1774  & 424  & 37.3  & 158 \tabularnewline
{\small{}Baalback} & 877  & 2146  & 94  & 40.6 & 38 \tabularnewline
{\small{}Hermel} & 86  & 305  & 42  & 47.1  & 20 \tabularnewline
\hline
\end{tabular} {\small\par}
\par\end{raggedright}

\end{table}
\par\end{flushleft}

\section {Conclusion}

 In this paper we introduced the Moran's $I$ index with its associated $z$-score and $p$-value to study the spatial autocorrelation of registered new infections of COVID-19 in Lebanon. We introduced six different cases of parameterization of the spread related to adjacency, proximity, population, population density, poverty rate and poverty density. We discovered that poverty rate is not statistically relevant to the spatial spread of the disease while geographic bordering, distance between district centers, number and density of residents and poverty density lead to clustering of the disease, with varying strengths and level of confidence since July and August through October. We also introduced methods to determine the geographic coordinates of the mean center of the infection, and determined this center since April 2020, and plotted its variations over time up until October.
The understanding of the spatial, demographic and geographic aspects of the disease spread over time provides an essential basis for the relevant authorities to take more efficient decisions of local and inter and intra-regional measures, thus contributing to increased social and health safety and security in the fight against the pandemic.

\section*{Appendix}

 The expected value of Moran's $I$ statistic is given by:
\[E[I]=\frac{-1}{N-1}\]

while its variance is defined as:

 \[V[I]=E[I^{2}]-E[I]^{2}\]
 where

\[E[I^{2}]=\frac{A-B}{\left(N-1\right)\left(N-2\right)\left(N-3\right)\left(\Sigma_{ij}W_{ij}\right)^{2}}\]

and $A$ and $B$ are given by:
\[
A=N\left[2\left(N^{2}-3N+3\right)\Sigma_{ij}W_{ij}^{2}-2N\Sigma_{i}\left(\Sigma_{j}W_{ij}\right)^{2}+3\left(\Sigma_{ij}W_{ij}\right)^{2}\right]
\]

\[
B=\frac{2\Sigma_{i}(X_{i}-\bar{X})^{4}}{\left(\Sigma_{i}(X_{i}-\bar{X})^{2}\right)^{2}}\left[\left(N^{2}-N\right)\Sigma_{ij}W_{ij}^{2}-2N\Sigma_{i}\left(\Sigma_{j}W_{ij}\right)^{2}+3\left(\Sigma_{ij}W_{ij}\right)^{2}\right]
\]
consequently, the $z_{I}$-score is given by $ z_{I}=\frac{I-E[I]}{\sqrt{V[I]}}$.

\section*{Acknowledgements}
The author has no acknowledgements.

\section*{Data Availability}
Raw data was obtained from publically available and cited resources. All used data and codes are available upon request.

\section{Introduction}

The spread of COVID-19 pandemic has practically affected the entire planet, and created enormous challenges on every aspect of human life and organization, starting with the health sector and with far reaching consequences on the economy, education, sports, transportation and politics. Since the first cases were registered in Wuhan, China in December 2019 \cite{WHO}, the global spatial dynamics of the infection have been changing as the disease swiftly moved towards the West \cite{Mathys1} into Europe then into the United States, South America, and eventually to the whole world, with nearly 38.1 million cases and 1.1 million deaths registered until October 12, 2020 \cite{Worldometer}.

Given the global geographic spread of the virus and the local wide spread in many countries, and the nature of the transmission of the virus, it is important to understand the spatial mechanisms of this spread and its dependence on proximity, demographics and social characteristics of infected areas. Spatial analysis provides a better understanding of the routes of transmission of infections \cite{Meng1}, consequently, it allows the decision-makers to draft and implement effective health and mitigation measures to reduce risks associated with the pandemic.

In Lebanon, the first case was registered on February 21, 2020 \cite{MOPH} and by October 12, $54624$ cases and $466$ deaths were registered \cite{DRM}. The first few weeks witnessed a relatively rapid increase but it sharply declined as a result of the strong mitigation measures enforced by the beginning of March. The lift of the international travel ban and the partial easing of measures led to the revival of higher spread rates since July. Only $1788$ cases were registered by July 1, 2020 before a sharp rise from July through October. The cases were mainly concentrated in Beirut, its suburbs and its neighboring areas in Mount Lebanon. On August 4, a huge explosion rattled the port of Beirut and destroyed thousands of houses and buildings in the surrounding areas. People were rushed into hospitals, with thousands of injuries recorded on that day \cite{BBC}. On such a horrible incident, hundreds of volunteers and civil defense teams were involved in rescue work for several days. The social distancing measures were largely neglected in such an emergency situation. The spread accelerated in the upcoming weeks, with sharp rise in Beirut and its surroundings and with a national widespread reaching all regions and major towns and cities \cite{Science}.

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.8]{Infection-spread}
\caption{A monthly map of the regional cumulative number of infections in Lebanese districts between April and October 2020.}\label{fig:1}
\end{center}
\end{figure}

\emph{Related Literature}: Spatial autocorrelation is the statistical analysis of data studied in space or in space-time aiming for the identification and estimation of spatial processes \cite{CliffOrd1, Anselin}. It has been implemented to study and analyze the spread of various diseases and infections including cancer, diabetes, SARS, influenza virus, COVID-19, etc... \cite{Glick, Hipp, SARS, PLOS, Kuwait}. Recent studies also inspected the effect of city size, population, transportation systems and demographics on the disease spread and its mortality rate \cite{PLOS1, Urban, Xie, Safety1, Safety2}.   The determination of the mean center of a population (centroid) was discussed in \cite{Aboufadel, Bachi, Barmore} and extending the concept to the determination of the mean center of wealth and infections allowed for a spatial analysis of the temporal dynamics of wealth distribution, economic growth and infectious diseases \cite{Mathys2}. The dynamics of the outbreak of COVID-19 in Lebanon and its reproduction number dynamics were studied in \cite{DeebJalloul1, DeebJalloul2, Mourad}.

In this paper, we study the clustering and spatial progression of new infections in Lebanon by applying the methods of spatial autocorrelation with different model parameterizations of geographic, demographic and social variables including  adjacency, proximity, population, population density, poverty rate and poverty density. Locating the mean center of the epidemic spread as a function of time is used to analyze the temporal geographic development of the spread. The obtained results provide a solid basis for the concerned policy makers to draw well-grounded and scientifically based local and regional measures that would contribute to controlling the infection spread.

 The paper is organized as follows: in section 2 we introduce the implemented analytic mathematical and statistical methods and tools. Results are presented and discussed in section 3, and section 4 concludes the paper.

\section{Analytic methods and tools}

\subsection{Moran's I index}

Moran's $I$ index is an inferential statistic used to measure the spatial autocorrelation based both on  locations and feature values simultaneously. It is defined as \cite{CliffOrd1}:
\begin{equation}
 I=\frac{N\Sigma_{ij}W_{ij}(X_{i}-\bar{X})(X_{j}-\bar{X})}{\Sigma_{ij}W_{ij}\Sigma_{i}(X_{i}-\bar{X})^{2}}
\end{equation}

where $W_{ij}$ represents different types of adjacency between
region $i$ and region $j$, corresponding to different models of
infectious spread. $N$ is the number of regions under consideration
and $X_{i}$ represents the number of new daily infections
in district $i$. $\bar{X}$ is the average number
of new daily infections per region, and it is given by $\bar{X}=\frac{\Sigma_{i}X_{i}}{N}$.
The numerical outcome of $I$ falls between $-1$ and $1$ and it indicates whether a distribution is dispersed, random or clustered. A value of $I$ close to $0$ indicates a random distribution, while positive values  indicate clustered spatial distribution and negative values indicate dispersion. Larger values of $|I|$ nearer to $1$ mean stronger clustering (positive $I$) or stronger dispersion (negative $I$).

The $z_{I}$-score associated to this statistic is defined by:

\begin{equation}
z_{I}=\frac{I-E[I]}{\sqrt{V[I]}}
\end{equation}

where the expected value $E[I]$ and the variance $V[I]$ are defined in the Appendix. The $z$-score or the corresponding $p$-value of the statistic are used to reject the null hypothesis and eliminate the possibility of a random pattern leading to the obtained value of the Moran $I$ statistic.

In this paper, we take a $95\%$ confidence level corresponding to $|z_{I}|>1.96$ or equivalently to $p<0.05$ in order to confirm the outcome of clustering or dispersion of our spatial data indicated by $I$. In this case we say that the $p$-value is statistically significant, and based on the value of $I$ we can determine the pattern of the distribution.

We consider a model with six different cases of parameterization of the adjacency matrix $W_{ij}$ corresponding to geographic adjacency (case I), proximity (case II), population (case III), population density (case IV), poverty rate (case V) and poverty density (case VI). Table \ref{table:1} summarizes relevant data from the Lebanese districts.

In  casel I, we take $W_{ij}=1$ for districts sharing common borders and $W_{ij}=0$ otherwise, while in case II we determine $W_{ij}=\frac{1}{d_{ij}}$ where $d_{ij}$ is the driving distance between the administrative centers of regions $i$ and $j$. Those two cases study the effect of administrative adjacency and the distance proximity of different districts on the geographic clustering of new infections in Lebanon.

In case III and case IV, we apply the methods used in (\cite{Meng1, Kang}) to analyze the effects of population and population density on the spread of the disease since the virus is carried by people and its spread is supposed to be related to their population and interaction. We sort the districts by the number of their residents (obtained from \cite{CAS})  and then by the density of their residents relative to their areas. In these two cases, districts of consecutive populations and population densities  are assigned a factor of $W_{ij}=1$, and $W_{ij}=0$ otherwise. This provides a statistic about the clustering of infections according to population and population density respectively.

Lastly, in cases V and VI, we introduce  new parameters, namely the poverty rate and the poverty density in different districts and we analyze their effect on infection clustering. We sort the distritcs by their rates of poverty  and poverty density  \cite{CAS} and assign  $W_{ij}=1$ for regions of consecutive order of poverty rate or poverty density, and $W_{ij}=0$ otherwise, in a similar methodology to cases III and IV in order to infer the effect of poverty rate and density on the geographical patterns of infection spread.

\subsection{Mean center of infection}

The mean center of infection (henceforth MCI) is a geographic location that represents the weighted mean of the positions of infected individuals on the surface of Earth, assumed to be spherical. Assigning the value of Earth's radius to unity, the two spherical coordinates that determine the unique position of a point are its latitude $\lambda_{i}$ and longitude $\phi_{i}$.  The latitude is a measurement of location north or south of the equator while the longitude is a measurement of location east or west of the prime meridian at Greenwich, UK.

The Cartesian position vector $\overrightarrow{r_{i}}=(x_{i},y_{i},z_{i})$
is related to spherical coordinates with unit radius by the relations \cite{Coordinates}:

\begin{equation}
\begin{cases}
x_{i}=\cos\lambda_{i}\cos\phi_{i}\\
y_{i}=\sin\lambda_{i}\cos\phi_{i}\\
z_{i}=\sin\phi_{i}
\end{cases}
\end{equation}

We denote the district number of infections (new or cumulative) by $X_{i}$ as defined above, and the Cartesian positions of the administrative centers by $(x_{i},y_{i},z_{i})$ . Then, the Cartesian position of the weighted mean of infections
$\overrightarrow{\hat{r_{i}}}$ is given by:

\begin{equation}
\begin{cases}
\hat{x}=\frac{\Sigma_{i}X_{i}x_{i}}{\Sigma_{i}X_{i}}\\
\hat{y}=\frac{\Sigma_{i}X_{i}y_{i}}{\Sigma_{i}X_{i}}\\
\hat{z}=\frac{\Sigma_{i}X_{i}z_{i}}{\Sigma_{i}X_{i}}
\end{cases}
\end{equation}

As suggested by \cite{Aboufadel}, the precise position on the surface of a sphere can
be determined from the normalized position vector defined by $\overrightarrow{\bar{r_{i}}}=(\bar{x},\bar{y},\bar{z})=\frac{\overrightarrow{\hat{r_{i}}}}{\mid\overrightarrow{\hat{r_{i}}}\mid},$leading
to

\begin{equation}
\begin{cases}
\bar{x}=\frac{\hat{x}}{\sqrt{\hat{x}^{2}+\hat{y}^{2}+\hat{z}^{2}}}\\
\bar{y}=\frac{\hat{y}}{\sqrt{\hat{x}^{2}+\hat{y}^{2}+\hat{z}^{2}}}\\
\bar{z}=\frac{\hat{z}}{\sqrt{\hat{x}^{2}+\hat{y}^{2}+\hat{z}^{2}}}
\end{cases}
\end{equation}

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.55]{Neighbourhood}
\caption{The figure shows Moran's I index and its corresponding p-value for cases I and II accounting for adjacency and proximity of new infections registered in Lebanese districts. The test shows that since July 2020, there is clear geographic clustering of infections. }\label{fig:2}
\end{center}
\end{figure}

\begin{figure}[h]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.55]{Demographics}
\caption{Moran's I index and its corresponding p-value for cases III and IV accounting for population and population density of different districts. Population density contributes for strong clustering since August 2020, while the districts' population effect is weaker. } \label{fig:3}
\end{center}
\end{figure}

Consequently, we can recover the spherical position of the mean center of infections by calculating the mean latitude and longitude as:

\begin{equation}
\begin{cases}
\bar{\phi}=\sin^{-1}\bar{z}\\
\bar{\lambda}=\tan^{-1}\frac{\bar{y}}{\bar{x}}
\end{cases}
\end{equation}

The latitude and the longitude can be located and plotted on maps and geographic information systems. We employ the spherical coordinates of geographic locations of the capitals of the $26$ administrative districts in Lebanon and the number of daily and cumulative infections in each region in order to determine the daily MCI accordingly. This provides a tool to analyze the temporal dynamics of the mean geographic spread of the disease.

\section{Results and discussions}

The determination of the Moran's $I$ index and its corresponding p-value for the effect of adjacency and proximity of cases I and II on the clustering of daily new infections of COVID-19 in Lebanon shown in Figure (\ref{fig:2}), leads to the conclusion that since July 2020, there is strong clustering of infections in regions sharing common borders and among nearer regions. There were only few days when new infections were not clustered in adjacent regions, and only one day where distance was not shown to be a detrimental effect in the spatial spread of new cases. The maximum value of Moran's $I$ reached $0.660$ for case I and $0.380$ for case II indicating a high level of geographic clustering of the disease spread since July. The infections before July had a high $p$-value, indicating a high probability for random geographic spread.

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.55]{Poverty}
\caption{Moran's I index and its corresponding p-value for cases V and VI accounting for regional poverty rate and poverty density. The poverty rate is not a decisive factor for spatial spread but poverty density contributes to spatial clustering since the end of August 2020. }\label{fig:4}
\end{center}
\end{figure}

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.55]{Centroid}
\caption{This figure is a plot of the latitude and longitude of the weighted geographic center of COVID-19  infections in Lebanon. The upper graphs represent the temporal progression of the cumulative number of infections while the lower graphs represent that of the new daily cases.}\label{fig:5}
\end{center}
\end{figure}

The results of the spatial spread dynamics in relation to population and population density adjacency as shown in Moran's $I$ and $p$-value of cases III and IV depicted in  Figure (\ref{fig:3}) reveal that the spread was not clustered with respect to the regional population until late August 2020, where it started achieving a positive value of $I$ with $p<0.05$ indicating spatial clustering between regions of adjacent population rank, with several days showing a probability of random spread. The maximum attained $I$ was $0.666$. However, the statistics for districts with adjacent rank of population density show very strong spatial clustering since the middle of August with $I$ attaining a maximum value of $0.832$, which is the highest among all six studied cases.

The results of case V  (Figure \ref{fig:4}) show that the spatial spread cannot be attributed to adjacent ranking of poverty rates among the districts since the p-values remain above the $5\%$ level of confidence up until October 2020, hence no spatial clustering occurs. But when we consider the poverty density in case VI, we obtain positive values for Moran's $I$ since the end of August, with $p<0.05$ except for five days. Hence, spatial clustering among regions with adjacent ranking of poverty density occurs. The maximum attained $I$ in this case is $0.666$.

In comparison, we find out that clustering of new infections occurs starting on different dates between July and August for all considered cases except for case V corresponding to district populations. The strongest level of spatial clustering (highest $I$) occurs for model IV of population density after mid-August, while clustering associated to geographic adjacency and proximity (cases I and II) has the longest time span (since early July) and the highest levels of confidence.

\begin{figure}[t]
\hspace*{-2cm}
\begin{center}
\includegraphics[scale=0.65]{MCI}
\caption{The figure shows the variation of the geographic location of the MCI on a monthly basis between April 1 and October 1, 2020, together with the mean center of the population.}\label{fig:6}
\end{center}
\end{figure}

The location of the MCI was determined as a function of time as shown in   Figure (\ref{fig:5}). The mean latitude and longitude of the infection were determined according to the methods described in equation (6). The location of the cumulative MCI is plotted on the geographic map of Lebanon during the same period in  Figure (\ref{fig:6}), together with the mean center of population of the country. It started near the city of Jounieh, north-west of the mean center of population, but it has moved southward since May through August, where it started moving northward again. The location of the MCI of new infections was quite geographically distributed before July as the lower plot of  Figure (\ref{fig:5}) shows, before becoming more homogenous afterwards.

The reproduction number $R$ and the rate of the infection spread correlate with people's mobility \cite{Safety3}. Geographic clustering occurs because people's motion and local travel is higher in their close neighborhoods, especially in a country like Lebanon where with the absence of national public transportation throughout the country \cite{Transports} diminishes nationwide mobility. Higher levels of social interaction among people in dense regions also contribute to the spread of the disease, and this has shown the strongest clustering effect. The study of the spread of the infection allows the relevant authorities to draw appropriate country-specific and regional measures to curb the spread.

\begin{flushleft}

\begin{table}[t]
\begin{raggedright}
{\small{}
\caption{The table shows the distribution of the cumulative number of cases among the 26 Lebanese districts on October 12, 2020, with their respective populations, population densities, poverty rates and poverty densities.} \label{table:1}
\smallskip
\begin{tabular}{lc lc lc lc lc lc}
\hline
{\small{}Region Name} & {\small{}Number of cases} & {\small{}Population } & {\small{}Population denisty} & {\small{}Poverty} & {\small{}Poverty density}\tabularnewline
{\small{}} & {\small{} (Cumulative)} & {\small{} $\times100$} & {\small{}(resident/$km^{2}$)} & {\small{}rate (\%)} & {\small{}(resident/$km^{2}$)}\tabularnewline
\hline
{\small{}Akkar} & 1171 & 3204 & 418  & 38.4 & 38 \tabularnewline
{\small{}Minieh-Denniyeh} & 723  & 1408  & 389  & 48.6  & 189 \tabularnewline
{\small{}Tripoli} & 4198 & 2438 & 9030 & 31.7 & 2862 \tabularnewline
{\small{}Zgharta} & 854 & 877 & 399  & 25  & 100 \tabularnewline
{\small{}Koura} & 556 & 846 & 489 & 14.3 & 70 \tabularnewline
{\small{}Bcharre} & 115 & 221  & 140  & 13.4  & 19 \tabularnewline
{\small{}Batroun} & 341 & 589 & 212 & 5.5 & 12 \tabularnewline
{\small{}Jbeil} & 954 & 1295 & 301  & 12.7  & 38 \tabularnewline
{\small{}Kesrwan} & 1978 & 2605  & 762 & 18.5 & 141 \tabularnewline
{\small{}Meten} & 6139  & 5110  & 1928  & 20.4  & 393 \tabularnewline
{\small{}Beirut} & 6443  & 3417  & 17258  & 25.5 & 4401 \tabularnewline
{\small{}Baabda} & 7277  & 5538  & 2855 & 26.8  & 765 \tabularnewline
{\small{}Aley} & 3047 & 3008 & 1144  & 29.4  & 336 \tabularnewline
{\small{}Chouf} & 1965  & 2770  & 560  & 24.2  & 135 \tabularnewline
{\small{}Jezzine} & 125 & 321  & 133  & 21.9  & 29 \tabularnewline
{\small{}Saida} & 2472  & 2966  & 1079  & 19.1  & 206 \tabularnewline
{\small{}Nabatieh} & 683  & 1802  & 593  & 28.2  & 167 \tabularnewline
{\small{}Sour} & 1023 &  2557 & 933  & 30.3 & 283 \tabularnewline
{\small{}Bent Jbeil} & 331 & 962  & 364  & 22.9  & 83 \tabularnewline
{\small{}Marjeyoun} & 214 & 740  & 279  & 24.2  & 68 \tabularnewline
{\small{}Hasbaya} & 76 & 287  & 108  & 23.9  & 26 \tabularnewline
{\small{}Rachaya} & 79 & 338  & 62  & 16  & 10 \tabularnewline
{\small{}West Beqaa} & 474  & 864 & 184  & 25.5  & 47 \tabularnewline
{\small{}Zahleh} & 2367 & 1774  & 424  & 37.3  & 158 \tabularnewline
{\small{}Baalback} & 877  & 2146  & 94  & 40.6 & 38 \tabularnewline
{\small{}Hermel} & 86  & 305  & 42  & 47.1  & 20 \tabularnewline
\hline
\end{tabular} {\small\par}
\par\end{raggedright}

\end{table}
\par\end{flushleft}

\section {Conclusion}

 In this paper we introduced the Moran's $I$ index with its associated $z$-score and $p$-value to study the spatial autocorrelation of registered new infections of COVID-19 in Lebanon. We introduced six different cases of parameterization of the spread related to adjacency, proximity, population, population density, poverty rate and poverty density. We discovered that poverty rate is not statistically relevant to the spatial spread of the disease while geographic bordering, distance between district centers, number and density of residents and poverty density lead to clustering of the disease, with varying strengths and level of confidence since July and August through October. We also introduced methods to determine the geographic coordinates of the mean center of the infection, and determined this center since April 2020, and plotted its variations over time up until October.
The understanding of the spatial, demographic and geographic aspects of the disease spread over time provides an essential basis for the relevant authorities to take more efficient decisions of local and inter and intra-regional measures, thus contributing to increased social and health safety and security in the fight against the pandemic.

\section*{Appendix}

 The expected value of Moran's $I$ statistic is given by:
\[E[I]=\frac{-1}{N-1}\]

while its variance is defined as:

 \[V[I]=E[I^{2}]-E[I]^{2}\]
 where

\[E[I^{2}]=\frac{A-B}{\left(N-1\right)\left(N-2\right)\left(N-3\right)\left(\Sigma_{ij}W_{ij}\right)^{2}}\]

and $A$ and $B$ are given by:
\[
A=N\left[2\left(N^{2}-3N+3\right)\Sigma_{ij}W_{ij}^{2}-2N\Sigma_{i}\left(\Sigma_{j}W_{ij}\right)^{2}+3\left(\Sigma_{ij}W_{ij}\right)^{2}\right]
\]

\[
B=\frac{2\Sigma_{i}(X_{i}-\bar{X})^{4}}{\left(\Sigma_{i}(X_{i}-\bar{X})^{2}\right)^{2}}\left[\left(N^{2}-N\right)\Sigma_{ij}W_{ij}^{2}-2N\Sigma_{i}\left(\Sigma_{j}W_{ij}\right)^{2}+3\left(\Sigma_{ij}W_{ij}\right)^{2}\right]
\]
consequently, the $z_{I}$-score is given by $ z_{I}=\frac{I-E[I]}{\sqrt{V[I]}}$.

\section*{Acknowledgements}
The author has no acknowledgements.

\section*{Data Availability}
Raw data was obtained from publically available and cited resources. All used data and codes are available upon request.