\section{Introduction}
Almost twenty years after the observational evidence of cosmic acceleration the cause of this phenomenon, dubbed
 "dark energy" (DE), remains an open issue which challenges the foundations of theoretical physics: whenever the DE is the simplest option - a cosmological constant $\Lambda$ or some alternative theory of gravity explaining cosmological dynamics and large scale structure
\cite{Capozziello:2002rd,Nojiri:2010wj,Tsujikawa:2010zza,Capozziello:2011et,Cai:2015emx,Nojiri:2017ncd,Carlesi:2016yas,Odintsov:2020nwm}. Modified theories of gravity are theoretically and observationally appealing, and ever-increasingly sensitive and precise upcoming surveys provide the exciting possibility of robustly testing them against observational data. Among these alternatives,  extensions of General Relativity can be considered a straightforward and natural approach to retain positive results of Einstein theory and eventually to extend it at infrared and ultraviolet scales (see \cite{Capozziello:2011et} for a detailed discussion).

Beside the cosmological systems giving the probes for DE, the LG of galaxies can be used to track DE effects  \cite{Hoffman:2007fu,Erdogdu:2009zb,Eingorn:2012dg,Eingorn:2012jm,Partridge:2013dsa,Gonzalez:2013pqa,McLeod:2016bjk,McLeod:2019cfg}. So far, many tries to test the cosmological constant effect in the LG dynamics have been considered. Here we extend the research for extended gravity models (such as $f(R)$ gravity) where Yukawa-like potentials naturally emerge in the low energy limit \cite{Stabile, Capozziello:2020dvd,Cardone:2011ze}.  {The only consistent $f(R)$ models which evade the solar system constraints have a chameleon mechanism gravity \cite{Brax:2008hh}.} Assuming the LG mass from different simulations, it is possible to  yield bounds for the Yukawa parameters. This research comes into a very large debate related to the dark matter.

In fact, dark matter is one of the most profound unsolved phenomena in modern astrophysics and cosmology. The standard approach describes dark matter as cold massive particles dubbed Weakly Interacting Massive Particles (WIMPs) \cite{Tao:1989xn,Morales:2002ud,Iocco:2012jt,Conrad:2014tla,Rott:2012gh,Baudis:2013eba,DRUKIER:2013lva,daSilva:2014qba,Cui:2015eba,Arcadi:2017kky,Queiroz:2017kxt}, axions \cite{Peccei:1978fx,Davier:1987dv,Murayama:1998jb,Kim:1998sy,Kim:1999ia,vanBibber:2001ud,Geralis:2009ue,Kim:2009xp,Pajer:2013fsa} or very light axion-like particles \cite{Masso:2002ip,Galanti:2019sya,Ertas:2020xcc}. The presence of dark matter in galaxies is observed from different measurements. The basic one is the mismatch between the predicted Keplerian velocity of orbiting stars in galaxies and the measured one \cite{Rubin:1980zd,Begeman:1991iy}. The measured velocity for large distances is approximately constant. In addition, this constant velocity is related to the luminous mass through the Tully-Fisher relation \cite{Zwaan:1995uu,McGaugh:2000sr,TorresFlores:2011uc,McGaugh:2011ac,Chen:2019ftv}.  Other systems where  dark matter is present are the galaxies. Despite of these relevant astrophysical evidences, no new fundamental particle  has been detected till now neither by direct  nor by indirect detection to address the dark matter issue. In this perspective, instead of searching for new matter ingredients, extragalactic and large scale structure dynamics could be addressed by extending gravitational sector (see e.g. \cite{Capozziello:2017rvz,Capozziello:2020dvd, Napolitano:2012fp}).

An important approach in this debate is constituted by the Modified Newtonian Dynamics (MOND) \cite{Milgrom:1983ca,Bekenstein:1984tv,Sanders:2002pf} which is an alternative explanation to the flat rotation curves for galaxies not implying any new matter ingredient. By violating Newton's second law at low accelerations, the Tully-Fisher relation is recovered, without introducing additional dark matter. The Tensor Vector Scalar theories (TeVeS) is a covariant theory, which in the low energy limit, produces MOND \cite{Bekenstein:2004ne}. Specifically, in \cite{Capozziello:2017rvz}, it is  shown that $f(R)$ gravity yields the flat rotation curves of galaxies recovering MOND in the weak field limit. In the second part of this paper, we investigate the LG dynamics, assuming that dark matter is an effect emerging from some $f(R)$ gravity adopting the same philosophy as in \cite{Capozziello:2017rvz} and we show that the LG gives prediction for the ratio between and dark matter and dark energy that could be tested with the observations.

The structure of the paper is the following: Section \ref{sec:form} formulates the timing argument (TA) and the method we use. Section \ref{sec:FRDE} discusses the effect of Yukawa potential on the LG mass and constraints the Yukawa parameters. Section \ref{sec:FRDM} is devoted to modified gravity the capable of mimicking dark matter in the LG system. Here, we calculate the predicted mass of the LG. Section \ref{sec:dis} summarizes the results.

\section{The total Mass from the Timing Argument}
\label{sec:form}
The Cosmological Constant  is considered to govern dynamics at cosmological scales \cite{Chernin:2000pq,Baryshev:2000kw,Chernin:2001nu,Karachentsev:2003eh,Chernin:2003qd,Teerikorpi:2005zh,Chernin:2009ms,Chernin:2006dy,Teerikorpi:2010zz,Chernin:2015nna,Silbergleit:2019oyx}. Hence we include the Cosmological Constant contribution in our analysis. The center of mass coordinate system is defined by the relative distance $r$ and the relative velocity $v$. The relative distance variation reads \cite{Emelyanov:2015ina,Carrera:2006im}:
\begin{equation}\label{ENL}
\ddot{r} = \frac{l^2}{r^3}-\frac{GM}{r^2} +  \frac{1}{3}\Lambda c^2 \, r,
\end{equation}
where $\vec{l} =  \vec{r} \times \vec{v}$ and $M$ is the total mass of MW and $M_{31}$. Fig. \ref{fig:potential} presents the corresponding effective potential with and without $\Lambda$. It is possible to see that $\Lambda$ is important even in these short scales \cite{Eingorn:2012dg,Eingorn:2012jm,Partridge:2013dsa,Gonzalez:2013pqa,McLeod:2016bjk,McLeod:2019cfg,Benisty:2019fzt,Lemos:2020vhj}. There are two extreme points: The closer one is a minimum and the other is a maximum. For large distances, the repulsive force arising from the Cosmological Constant makes impossible any bound structure so the Cosmological Constant gives us a constraint on the maximum size of the bound systems.

\begin{figure}[t!]
 	\centering
\includegraphics[width=0.49\textwidth]{pot.pdf}
\caption{\it{ {The effective potential vs. radius for two-body system, assuming $M = 4\cdot 10^{12} \, M_{\odot}$, for different cases: \textbf{dashed - blue:} Newtonian potential, \textbf{orange:} Newtonian potential with Cosmological Constant (from Planck 2018 value), \textbf{green:} Yukawa type potential (with the values $\alpha = 0.1, m_g = 0.01_{Mpc^{-1}}$) and with the Cosmological Constant.}}}
 	\label{fig:potential}
  \end{figure}
\begin{figure}
 	\centering
\includegraphics[width=0.47\textwidth]{SpectFullGR.pdf}
\caption{\it{ {The evaluated distance between Milky Way and $M\text{31}$ at $t=0$ for different TA masses of the LG. The model predicts the correct mass when $r(t=0) = 0$. The first minimum ($4-5\, 10^{12} M_{\odot}$) points the predicted mass according to the TA. The second minimum points the predicted mass with $\Lambda$. }}}
 	\label{fig:massdistance}
  \end{figure}

The Timing Argument (TA) \cite{1989ApJ345108P,1981Obs101111L} for the Local Group is one of the historical probes for the missing mass problem. The TA assumes that Milky Way ($MW$) and Andromeda ($M31$) have been approaching each other despite cosmic expansion. In its simplest version, the LG roughly consists of the $MW$ and $M31$ as two isolated point masses. Briefly after their formation (for the sake of simplicity the "big bang"),  these two galaxies must have been in the same place with zero separation. Due to the Hubble expansion, these two galaxies moved apart. After a couple of billion years, they slowed down and then moved towards each other again, as a consequence of the gravitational pull. We restrict ourselves to a simple basis of the TA as an isolated, two dimensional system. Galaxies are modeled as point masses. Moreover, as the galaxy pairs are isolated, there are no external gravitational fields. The "initial condition" of the model is $r = 0$, which refers to the "big bang" condition. In order to calculate the mass of the LG, we change the direction of time by considering the opposite direction of the measured velocity now. We evaluate the measured distance of $M31$ to obtain what should be the distance at the "big bang" for different LG masses.

Fig.~\ref{fig:massdistance} presents the distance for the big bang for different masses of LG. When the curve approaches minimum ($r \rightarrow 0$) yielding the predicted mass. The minimum recognizes the mass of the LG as the model predicts. We use the units which can get the current values of the distance and the velocities of $\text{M31}$:
\begin{equation}
\begin{split}
r_{m_{31}} = 0.77 \pm 0.04 \, {\text{Mpc}}, \\ v_{rad} = -109.3 \pm 4.4 \, {km/sec},  \\
v_{tan}(t_u) = 82^{+38}_{-35} \, {km/sec},
\end{split}
\end{equation}
  {based on \cite{vanderMarel:2012xp,2021MNRAS.507.2592S}. Here $v_{rad}$ is the radial velocity towards the $MW$. The cosmological constant is determined by the latest Planck measurements \cite{Aghanim:2018eyx}. When we estimate the initial distance at the big bang, we use the age of the Universe. We use the initial conditions as a Gaussian prior. }

 {Given the simplicity of the Timing Argument and the assumption that the MW and M31 are point masses with constant masses, we have to include some corrections to get unbiased estimates of the true LG mass. This correction arises from other LG members and from the halos around the MW and M31 galaxies \cite{Gupta:2012rh,Bregman:2018dne,Bogdan:2022pbc}. We correct the inferred LG masses with the cosmic bias effect based on \cite{Hartl:2021aio,Benisty:2022ive} with a Gaussian prior of $0.62 \pm 0.2$. The mass of the LG without the Cosmological Constant is: $M_{LG} = \left( 3.28 \pm 0.31\right) \cdot 10^{12} \, M_{\odot},$ and with the Cosmological Constant presence we get: $M_{LG} = \left(3.86 \pm 0.45\right) \cdot 10^{12}  \, M_{\odot}.$ The posterior distributions of the LG mass in presented in Fig~\ref{fig:massDis}. It is possible to see the correlation between the LG mass and the distance of M31 towards us. Therefore for Yukawa type potential that depends on $r$ the predicted mass would be changed and could be constraint based on the history of the LG. }

\begin{figure}[t!]
 	\centering
\includegraphics[width=0.39\textwidth]{massHisGRLam.pdf}
\caption{\it{The predicted masses from the Timing Argument + Cosmic Bias considering the absence and the  presence the Cosmological Constant. The total mass distribution has a correlation to the distance of $M31$ towards us. Therefore for Yukawa potential that depends  }}
 	\label{fig:massDis}
  \end{figure}

\section{Yukawa Potential Constraints}
\label{sec:FRDE}
One of the simplest way to test gravity is to constrain the existence of a Yukawa additional term, such that the one particle Newtonian potential assumes the form:
\begin{equation}
\Phi(r) = \frac{G M}{r}\frac{1 + \alpha \, \exp\left[ - m_{grav} r\right]}{1+\alpha},
\end{equation}
where $M$ is the total mass of the system, $G $ is the Newtonian constant and $r$ is the separation between the two objects. For the Yukawa couplings: $\alpha$ is the Yukawa strength and $m$ is the Yukawa mass. For $\alpha$ goes to zero the Yukawa interaction is reduced into the Newtonian one. For $m_{grav}$ goes to zero the Newtonian constant is modified to $G \left(1 + \alpha\right)$. Such a correction arises whenever the force is mediated by a scalar particle of mass $m$. Therefore testing gravity is equivalent to testing the existence of a fifth force of scalar nature.

Fig~\ref{fig:potential} shows the potential for different values of $\alpha$ and $m_{grav}$. The contribution of the Cosmological Constant effects further the scale of the LG (more then $\sim 1 Mpc$). However, the perturbations of $f(R)$ affect  scales closer than the LG scale. This model cannot replace the Cosmological Constant effect and the latter needs to be added as an additional perturbation.

\begin{figure}[t!]
 	\centering
\includegraphics[width=0.49\textwidth]{DiffAlpha.pdf}
\caption{{\it{The dependence of the LG mass vs. different values of $\alpha$ and the $m_{\text{grav}}$ in $Mpc^{-1}$. Negative values of $\alpha$ increasing the total mass while positive values of $\alpha$ reducing the total mass. The colors represent different values of $\alpha$: black - $\alpha = 0$, blue - $\alpha = \pm 0.01$, red - $\alpha = \pm 0.1$ and green - $\alpha = \pm 0.5$.}}
\label{fig:alphaMass}
	\end{figure}

As shown in several studies, Yukawa-like corrections emerges as common  features for several metric theories of gravity \cite{Stabile, Mairi}. In particular for theories containing higher-order curvature invariants or scalar-tensor theories \cite{Capozziello:2010zz}. The field equations for $f(R)$ gravity reads:
\begin{equation}\label{fe}
f'(R)R_{\mu\nu}-\frac{1}{2}f(R)g_{\mu\nu}-f'(R)_{;\mu\nu}+g_{\mu\nu}\Box
f'(R)= T_{\mu\nu}\,,
\end{equation}
with the  trace
\begin{equation}\label{fetr}
3\Box f'(R)+f'(R)R-2f(R)= T \,.
\end{equation}

where $\kappa$ is the gravitational coupling and $T_{\mu\nu}$ is the matter energy-momentum tensor. We are interested in external solutions so we can ignore the matter contribution.  In the weak field limit, we can perturb the metric tensor with respect to the Minkowski background:
\begin{equation}
g_{\mu\nu}\,=\,\eta_{\mu\nu}+h_{\mu\nu},
\end{equation}
where $|\eta_{\mu\nu}|\ll |h_{\mu\nu}|$. For an analytic
 $f(R)$ Lagrangian expandable in Taylor series
\begin{eqnarray}\label{sertay}
\begin{aligned}
f(R)=\sum_{n}\frac{f^n(R_0)}{n!}(R-R_0)^n\simeq\\
f_0+f'_0R+f''_0R^2+f'''_0R^3+...
\end{aligned}
\end{eqnarray}
where the prime indicates derivatives with respect to $R$,
the field equations, in the post-Newtonian limit, up to ${\mathcal O}(4)$ order. By adopting the spherical symmetry, the metric can be recast \cite{Capozziello:2012ie,DeMartino:2018yqf,DeLaurentis:2018ahr,DeLaurentis:2018udw,Cardone:2011ze,Napolitano:2012fp,Capozziello:2008ny,Hees:2017aal,Zakharov:2018cbj,Capozziello:2020dvd}:
\begin{equation}
\label{Edd}
    g_{tt}=1-\frac{\Phi(r)}{c^2}\,,\;\;\;\;\;\;  g_{rr}=1+\frac{\Psi(r)}{c^2}\,.
\end{equation}
 {where the functions $g_{tt}$ and $g_{rr}$ give two gravitational potentials:}
\begin{eqnarray}
\label{eq:PHI}\Phi(r) &=& -\frac{G M}{r} \frac{1+\alpha \,   e^{-m_{grav} r}}{1 + \alpha}, \\
\label{eq:PSI}
\Psi(r) &=& -\frac{G M }{r}\frac{1-\alpha \left(1 + m_{grav} r \right) e^{-m_{grav} r}}{1 + \alpha }\,.
\end{eqnarray}
 {The Yukawa strength and scale length read:}
\begin{equation}
\alpha = f_{0}' - 1, \quad m_{grav}^2\,=\,-f'_0/6f''_0,
\end{equation}
 {with a direct relation to the $f\left(R\right)$ function. The potentials coincide for $\alpha \rightarrow 0$ when Newtonian potential is recovered. }

In order to test the effect on the LG mass, we use Eq. ~(\ref{ENL}), where the Yukawa potential is replaces the Newtonian potential. Fig~\ref{fig:alphaMass} shows the numerical analysis and the predicted mass for different values of the LG. Negative values of $\alpha$ yields larger mass and positive values of $\alpha$ gives lower values of the mass.

\begin{figure}[t!]
 	\centering
\includegraphics[width=0.48\textwidth]{finalPost.pdf}
\caption{\it{The posterior distirubtion for the Yukawa parameters from the Local Group dynamics.   }	}
\label{fig:postFinal}
\end{figure}

In order to constraint the additional parameters we use the predicted mass from the latest LG mass prediction \cite{Benisty:2022ive}. We define the $\chi^2$ difference between the predicted mass and the simulated mass:
\begin{equation}
\chi^2 = \left(\frac{M_{GR + \Lambda} - M(\theta)}{\Delta M_{GR + \Lambda}} \right)^2,
\end{equation}
where $M(\theta)$ is the predicted mass from the set of parameters $\theta = \{r_0, v_r^{0}, v_t^{0}, \Lambda, \alpha, m_{grav} \}$ and $M_{GR + \Lambda}$ is the deduced mass from \cite{Benisty:2022ive}. Regarding the problem of $\chi^2$ minimization, we use nested Markov Chain Monte Carlo sampler as it is implemented within the open-source packaged $\text{Polychord}$ \cite{Handley:2015fda,Lewis:2019xzd}. Fig~\ref{fig:postFinal} shows the posterior distribution for the Yukawa potential. For a model with $m_{grav} = 0$ the largest deviation on $\alpha$ is:
\begin{equation}
\alpha = \left(2.807 +\- 3.01\right) \cdot 10^{-1},
\end{equation}
 {or: $\alpha < 0.58$. Assuming the Cassini satellite bound on $\alpha < 2 \cdot 10^{-5}$ with the LG likelihood could gives the constraint on the Yukawa mass:}
\begin{equation}
\begin{split}
 m_{grav} \preceq 5.094 \pm 2.855 \, Kpc^{-1} \\ = \left(3.265 \pm  1.830 \right) 10^{-26} \, eV/c^2.
\end{split}
\end{equation}
 {The error bars for the parameters are quite large but still indicate a difference from General Relativity. There is some degeneracy between the $\alpha$ and the $\lambda$ parameters that could resolve by assuming a prior on one parameter to get the other. The Yukawa mass is widely discussed in the literature \cite{deRham:2016nuf}. Since the LG scale is about $1 Mpc$. Therefore, it is clear why even the LG system that does not give a strong bound on the $\alpha$ parameter, still gives lower Yukawa mass, that corresponds to the LG large scale. }

 {The study of LG mass could be statistically stronger with another combined systems or with a complete simulation. In order to test the preference of the models we use the Bayesian Evidence (BE). The difference between the BE yields $ \Delta B_{ij} = 0.073 $ which indicates an \textit{indistinguishable preference} for General Relativity.}

\section{Modified Gravity or dark matter}
\label{sec:FRDM}
Let us discuss now whenever $f(R)$ gravity can replace the dark matter in galaxy clusters and in  LG in particular. Ref. \cite{Capozziello:2017rvz} addresses the dark matter in galaxies through the $f(R)$ model of a power\,-\,law case:
\begin{equation}
f(R)  = f_0 R^n\label{eq: frn}
\end{equation}
with $f_0$ a dimensional constant and $n$ any real number.
The choice of  power-law action \eqref{eq: frn} could appear non-natural in order to discuss deviations with respect to General Relativity.  If we assume small deviation with respect to General Relativity, and assume $n=1+\epsilon$ with $|\epsilon| \ll 1$, it is possible to re-write a first-order Taylor expansion as:
 \begin{eqnarray}\label{LOG}
R^{1+\epsilon}&\simeq & R+\epsilon R \,  {\rm log}R +O (\epsilon^2)\,,
 \end{eqnarray}
 and consider small deviations with respect to General Relativity.
The field equations can be set in the Einstein-like form \cite{Capozziello:2011et} as:
\begin{eqnarray}
R_{\mu \nu} - \frac{1}{2} g_{\mu \nu} R& = & \displaystyle{\frac{1}{f'(R)}}
\displaystyle{\Bigg \{ \frac{1}{2} g_{\mu \nu} \left [ f(R) - R
f'(R) \right ] + f'(R)_{; \mu \nu}} \nonumber \\ ~ & - &
\displaystyle{g_{\mu \nu} \Box{f'(R)} \Bigg \}} +
\displaystyle{\frac{T^{(m)}_{\mu \nu}}{f'(R)}} \label{eq:f-var2}\,.
\end{eqnarray}
 Primes denote derivative with respect to $R$. The  terms ${f'(R)}_{; \mu \nu}$ and $\Box{f'(R)}$ give rise to fourth order derivatives in $g_{\mu \nu}$. Taking into account the gravitational field generated by a point like source and solving the field equations (\ref{eq:f-var2}) in the vacuum case, we can write the metric as\,:
\begin{equation}
ds^2 = A(r) dt^2 - B(r) dr^2 - r^2 d\Omega^2 \label{eq: schwartz}
\end{equation}
where $d\Omega^2 = d\theta^2 + \sin^2{\theta} d\varphi^2$ is the spherical line element. Combining the $00$\,-\,vacuum component and the trace of the field equations (\ref{eq:f-var2}) in absence of matter, we get the equation\,:
\begin{equation}
\label{master-low} f'(R) \left ( 3 \frac{R_{00}}{g_{00}} - R
\right ) + \frac{1}{2} f(R) - 3 \frac{f'(R)_{; 0 0}}{g_{00}} = 0 \
.
\end{equation}
Using
Eq.(\ref{eq: frn}), Eq.(\ref{master-low}) reduces to\,:

\begin{equation}
\label{master-pla} R_{00}(r) = \frac{2 n - 1}{6 n} \ A(r) R(r) -
\frac{n - 1}{2 B(r)} \frac{dA(r)}{dr} \frac{d\ln{R(r)}}{dr} \ ,
\end{equation}
and the trace equation reads\,:

\begin{equation}
\Box{R^{n - 1}(r)} = \frac{2 - n}{3 n} R^n(r) \ . \label{eq:
tracebis}
\end{equation}
For $\epsilon = 0$, Eq.(\ref{eq: tracebis}) reduces to $R = 0$, which, inserted into Eq.(\ref{master-pla}), gives $R_{00} = 0$ and then the  standard Schwarzschild solution is recovered. Expressing $R_{00}$ and $R$ in terms of the metric (\ref{eq: schwartz}), Eqs.(\ref{master-pla}) and (\ref{eq: tracebis}) become a system of differential equations for  $A(r)$ and $B(r)$. A physically motivated
hypothesis is assuming:
\begin{equation}
A(r) = \frac{1}{B(r)} = 1 + \frac{2 \Phi(r)}{c^2} \label{eq:
avsphi}
\end{equation}
with $\Phi(r)$ the gravitational potential generated by a
point-like mass $m$ at the distance $r$ \cite{Capozziello:2006dp,Zakharov:2006uq}:
\begin{equation}
\Phi(r) = - \frac{G m}{2 r} \left [ 1 + \left ( \frac{r}{r_c}
\right )^{\beta} \right ] \label{equ02}
\end{equation}
\noindent where $r_c$ is an arbitrary parameter, depending on the
typical scale of the considered system and $\beta$ is a universal
parameter:
\begingroup
\setlength{\abovedisplayskip}{0pt}
\setlength{\belowdisplayskip}{0pt}
\begin{equation}
\beta = \dfrac{12 n^2 - 7n - 1 - \sqrt{36 n^4 + 12 n^3 - 83 n^2 + 50n
+ 1} }{6 n^2 - 4n + 2}.
\label{equ03}
\end{equation}
\endgroup
Clearly the gravitational potential deviates from the
Newtonian one because of the presence of the second term on the right hand side. The critical radius reads:
\begin{equation}
r_c = \sqrt{\frac{G M}{a_0}},
\end{equation}
where $a_0$ is the critical acceleration from MOND $a_0 = 10^{-10} m/s$. The parameter $r_c$ is an integration constant depending on $\beta$. Ref.~\cite{Capozziello:2007wc} shows that it is related to the fact that any $f(R)$ power-law  model is selected by  the existence of Noether symmetries with related conserved quantities. This fact can be considered a physical criterion to retain such models.

\begin{figure}[t!]
 	\centering
\includegraphics[width=0.51\textwidth]{SpectFullFRDM.pdf}
\caption{\it{The predicted Timing Argument mass of the LG for different values of $\beta$ with and without the Cosmological Constant. When $\beta \sim 0 $ we recover the Newtonian case. For some value of $\beta$ a minimal possible mass is obtained. }}
\label{fig:SpectFullFRDM}
	\end{figure}
\begin{figure}[t!]
 	\centering
\includegraphics[width=0.47\textwidth]{fraction.pdf}
\caption{\it{The statistical ratio between the the dark matter and the baryonic matter for the LG from the MONDian like potential prediction.}}
\label{fig:motionb}
	\end{figure}
 {It is possible to demonstrate that the potential (\ref{equ02}) can fit  the flat rotation curves in galaxies  like in MOND case \cite{Capozziello:2006dp}. As shown in \cite{Tula1,Tula2}, MOND dynamics can be fully recovered from $f(R)$ gravity. However here we adopt the TA for the potential (\ref{equ02}) for a binary system of galaxies starting from  Eq. (\ref{ENL}). Figure \ref{fig:SpectFullFRDM} shows the predicted mass of LG vs. the parameter $\beta$ that runs from zero to one, with or without a $\Lambda$ term and without the Cosmic Bias factor. The curve shows that there is a minimum predicted mass in a certain value of: }
\begin{equation}
\beta = 0.581^{+0.053}_{-0.081}.
\end{equation}
For any other value there is a duality for the mass with two different $\beta$ values.  {In order to test the distribution of the predicted mass, we take uniform prior on the parameter $\beta \in [0.;1.]$. The mass of  LG is $\left(2.35 \pm 0.21\right)\cdot 10^{12} \, M_{\odot}$ without $\Lambda$, and $ \left( 2.5 \pm 0.20  \right) \cdot 10^{12} \, M_{\odot}$ with $\Lambda$. }

The MONDian results that we consider here modify the potential while, in the original MOND  model, the formulation modifies the acceleration dependence. Refs.~\cite{McLeod:2019cfg,Benisty:2019fzt,Benisty:2020kys} test different formulations of MOND as a modified acceleration term and find that the prediction for  LG is around ten times lower. While for these formulations, MOND requires additional Past Encounters (PE) to get a reasonable number for the baryonic matter, the $f(R)$ version gives a case closer to the Newtonian one.

The ratio between the masses requires  a different analysis for MOND. In fact, MOND is a formulation that replaces dark matter. So the MONDian TA gives a prediction that all the mass is the baryonic matter alone. The Newtonian case predicts that the baryonic and an additional amount of dark matter constitute the total mass. The dark matter mass is calculated by the difference between the Newtonian prediction and the MONDian prediction. Figure
\ref{fig:motionb} shows the distribution for the ratio between the dark matter and the baryonic matter. The distribution yields the ratio:
\begin{equation}
\frac{Dark Matter}{Baryons} = 0.41^{+0.15}_{-0.38}
\end{equation}
with $1\,\sigma$ error.  {The ratio between dark matter and baryonic matter in our universe is around six and for the LG is about ten. The fact we get different number for the ratio shows that this function applies for one galaxy, but for two body system there is should be another function. This shows that the TA could test MONDian models of dark matter.}

\section{Discussion and Conclusions}
\label{sec:dis}
In this paper, we tested the LG of galaxies modeled as a two body problem  sensitive to the cosmological constant ($\Lambda$) contribution considered as an alternative   gravity effect. In the first part, we perturbed the dark energy effect assuming a Yukawa-like correction in the Newtonian potential. Specifically, we derived  the LG mass  from simulations and we constrained  the Yukawa strength. Still we get some ambiguity in the posterior distribution of Yukawa strength vs. the Yukawa length, but we showed that the mass is directly affected from the modified potential.

 {Constraints on $f(R)$ gravity are widely constraint in the literature \cite{Brax:2008hh}, from $\log_{10} f'(R_0) < -4.79$ in galaxy clusters \cite{Cataneo:2014kaa}, $\log_{10} f'(R_0) < -3$ from  GW 170817 \cite{Jana:2018djs} and $f'(R_0) < 3.7 \cdot 10^{-6}$ from the CMB \cite{Boubekeur:2014uaa}. Our constraint is not strong as these constraints ($f'(R_0) < 5.81 \cdot 10^{--1}$), however gives the constraint on $f(R)$ gravity for the first time from the LG dynamics.}

In the second part, we investigate the possibility that a minimal correction to General Relativity, i.e. $f(R) \sim R^{1+\epsilon}$ can replace the effect of dark matter as a MOND-like model. We find that there is a value of the parameter $\beta $, derived by $\epsilon$,  having a minimal value for the mass of  LG. Moreover, this particular potential predicts that the baryonic matter and the dark matter have, more or less, the same order of magnitude. A claim that could be tested with full a N-body simulation. As a conclusion,  we find that the LG dynamics can be considered a good test system for alternative theories of gravity which could be further improved with more realistic simulations and reliable data sets.

In general, modified gravity models of DM usually have some problems in fitting the CMB and other data sets and this is applies to $f(R)$ which mostly can fit galaxy rotation curves. Also in our case it is possible to see that one model of $f(R)$ gives a ratio of $0.41$ for the dark matter and baryonic matter ratio, and in order to fit for the experimental ratio the function possibly needs to be generelized.

\acknowledgments
D.B. thanks the Rothschild and the Blavatnik Fellowships for generous supports. This work has been supported by the European COST actions CA18108 and CA21136 and the research grants KP-06-N58/5. S.C. acknowledges the support of {\it Istituto Nazionale di Fisica Nucleare (INFN) iniziative specifiche QGSKY} and {\it MOONLIGHT2}.


\section{Main}\label{sec2}

A recent, high-impact development in the field of photonics has been the revelation of parity-time symmetric systems endowed with rich novel physics\cite{bender1999pt, bender1998real, el2018non, ruter2010observation, regensburger2012parity, jing2014pt,lin2011unidirectional,zhu2014p,benisty2011implementation,ding2015coalescence}. At the barycenter of the literature on PT-symmetry in photonics lie peculiar spectral degeneracies called the exceptional points\cite{peng2014parity,peng2014loss, cerjan2016exceptional,feng2013experimental,jing2017high,ding2016emergence,kang2016chiral,zhen2015spawning}. Originating from the realm of quantum mechanics, exceptional points (EPs) signify singularities at which the eigenvalues and eigenfunctions of non-Hermitian matrices cannot be described analytically\cite{heiss2012physics,kato1966analytic}. At the EP, two or more eigenvalues coalesce, simultaneously with the coalescence of associated eigenfunctions\cite{heiss1999phases,moiseyev2011non,el2007theory,pick2017general}. This theoretical construct has led to immense physical consequences in real experimental systems\cite{miri2019exceptional}. Previously discussed in the quantum theory of atomic and molecular resonances\cite{moiseyev1998quantum}, EPs are now demonstrated in carefully designed photonic systems motivated by the upgradation of systemic functionalities, such as, for instance, enhancement in output intensity\cite{feng2014single,hodaei2014parity,takata2021observing} and increased sensing capability\cite{hodaei2017enhanced,chen2017exceptional,hokmabadi2019non}. The ability of precise nanofabrication and on-demand gain/loss has allowed for the realization of several photonic platforms for the observation of EPs\cite{ozdemir2019parity,brandstetter2014reversing}. However, to our knowledge, there are no experimental reports of observations of EPs in disordered mesoscopic systems despite several investigations on non-Hermiticity in such systems\cite{vazquez2014gain,bachelard2022coalescence,davy2019probing,huang2021wave,weidemann2020nonhermitian,balasubrahmaniyam2020necklace,sahoo2022anomalous}.

The procedure of demonstration of exceptional points in photonic structures has typically followed a common theme. The gain or loss parameter in one cavity is systematically tuned while the emission spectrum and the spatial intensity in the sample is monitored. For example, in a coupled microring system, the authors fabricated two physically identical microrings at a pre-defined separation and coupling\cite{hodaei2014parity, hodaei2015parity}.  Optical excitation provided gain, and the loss in one of the cavities was continuously tuned by filtering out the pump energy. At the exceptional point, the lasing exclusively occurred over one microring, and the emission intensity was considerably higher than that from an individual microring pumped at par. This strategy has been adopted in various photonic systems \cite{takata2021observing, kim2016direct, chang2014parity, peng2014loss} that are carefully designed using physically separate cavities with tunable gain/loss. However, such a strategy is not feasible in an Anderson-localizing structure\cite{anderson1958absence,john1987strong,segev2013anderson,chabanov2000statistical,wiersma1997localization,sapienza2010cavity,schwartz2007transport,lahini2008anderson,riboli2011anderson} for several reasons.  Anderson-localizing systems are nondeterministic, wherein the occurrence of a resonance is probabilistic. In a situation of coupled  modes, both the cavities reside in the same physical structure and often involve common scatterers for distributed feedback. Naturally, calibrated introduction of gain/loss into individual localized modes is inconceivable. Hence, the conventional strategy of observing exceptional points is implausible in Anderson-localizing systems. In such a nondeterministic scenario, another route towards an EP can be adopted. The formation of an EP within a localizing system is a probabilistic event. Therefore, if a sufficiently large number of configurations can be generated and monitored, then an exceptional point can be expected with the right probability. Although easily said, this is a challenging prospect because the three ingredients , namely, (i) achieving Anderson localization, (ii) achieving coupled modes thereof, and (iii) addition of adequate gain/loss into the system, are all formidable experimental challenges. Nonetheless, this is the route we adopt in this work, and demonstrate Anderson-localization based lasing over exceptional points.

We use a technique which allows us to create and monitor thousands of disorder configurations\cite{joshi2020reduction, joshi2022anomalous}. The system is one-dimensional, and hence has an increased probability of Anderson localization and, therefore, coupled Anderson-localized modes.  Optical sources are distributed throughout the structure so that any localized modes existing at any spatial location are readily excited. We employ a simultaneous spectral, spatial and temporal diagnostic setup to completely characterize the emission from the microcavity array. The introduction of temporal diagnostics is a vital step towards identification of coupled modes. Usually, coupling is identified on the basis of spectral splitting. However, the spectra of Anderson localizing systems are random and multimodal in character. One cannot unambiguously distinguish a coupled-mode splitting from a pair of individual isolated modes with comparable spectral separation. On the other hand, the temporal signatures of the two cases are fundamentally different and provide conclusive evidence of splitting. Since temporal diagnostics are hitherto not reported at EP in lasing systems, we first simulate a deterministic prototype lasing structure and introduce the systematic temporal, spectral and spatial behavior at, and during the approach to, an EP. The EP-approach is characterized by quantum echoes (QE) that signify energy transfer between the coupled resonances, equivalent to Rabi oscillations in two-level systems\cite{dietz2007rabi,dembowski2004first}. The echoes vanish at the EP. Subsequently, we illustrate the experimental results in our system wherein the spectra exhibit coupled modes via spectral splitting and the temporal diagnostics simultaneously exhibit quantum echoes. The spatial intensity distributions exhibit the expected supermodes of the coupled system. The approach to an EP is evident in reduced splittings and lower oscillation frequencies, accompanied by increased output intensity. Eventually, emission peaks with largely enhanced intensity are shown, wherein the spatial intensity resides exclusively in one Anderson-localizing cavity. In these cases, the spectral profiles are seen to follow a square-Lorentzian profile, identifying an EP.

\begin{figure}[h!]
\centering
\includegraphics[width=1\textwidth]{simulation_arxiv.pdf}
\caption{\textbf{Quantum echo diagnosis of an exceptional point in a simulated, deterministic, model system}: (a) Temporal evolution of emission from two coupled microring resonators, as computed from finite-difference time-domain computations. Top panel: Situation for  spectral splitting of $\delta\lambda  = 3.6~$nm. Blue curve shows the echo in the emission at one wavelength, showing a period of $T$ ($=\frac{1}{\Omega}$, where $\Omega$ is the echo frequency). The echo in the emission at the other wavelength (orange curve) is out of phase, quantified by $\Delta T/T$. Middle panel: Dynamics for $\delta\lambda  = 1~$nm, where $T$ is increased, while the $\Delta T/T$ is reduced. Bottom panel: Behavior at the exceptional point (EP), i.e, $\delta\lambda  = 0~$nm, where there is no echo or oscillation. A single temporal profile is seen because of the eigenvalue coalescence. (b) Spectro-spatial observations corresponding to (a). Insets show intensity distribution in the coupled system averaged over the entire echo. The energy resides in both the resonators at large splittings, while it exists exclusively in the amplifying resonator at the EP. (c) Phase difference ($\Delta T/T$) in the two echoes (blue markers, left Y axis) reduces systematically with $\delta\lambda$, and is zero at the EP. Peak emission intensity (red markers, right Y-axis) $I_p$ diverges at the EP. Inset shows $\Omega$ with $\delta\lambda$. (d) Correspondence between the quantum echo behavior of (a) and PT-symmetry breaking and EP-formation. Circle markers mark the situations in (a) and (b).}
\label{fig:simring}
\end{figure}

We first illustrate the  quantum echo in a known EP system, namely, a coupled microring cavity. We simulate identical resonators with a diameter of $d = 2.0~\mu$m, thickness $t = 0.2~\mu$m and a separation $s = 0.4~\mu m$. The system was simulated using a finite-difference time-domain scheme (Lumerical) in the time domain, followed by full-wave finite-element computation (COMSOL) in the frequency domain. Optical gain or loss is introduced in the complex refractive index $n_{im}$ (see Methods for details). The $n_{im}$ in the second cavity is employed as the tunable parameter to approach an exceptional point, and is set to -0.01 (amplification) in the first cavity.  Fig\ref{fig:simring}(a) depicts the temporal behavior of the system, while (b) describes the corresponding eigenvalue behavior, along with the intensity distribution in the coupled resonators. The top panels of (a) and (b) depict the situation for identical gain, $n_{im} = -0.01$, in the two cavities. The resulting strong coupling realizes large splitting corresponding to a wavelength splitting $\delta\lambda$ of 3.6~nm. The emission clearly shows QE oscillations (top panel of Fig\ref{fig:simring}(a)) in the temporal profile. The blue and orange curves represent emission from the two cavities. The high QE frequency ($\Omega$) allows for multiple oscillations during the emission cycle, with out-of-phase oscillations from the two cavities. The corresponding split eigenvalues are shown in the top panel of (b), along with the intensity distribution (inset) in the cavities. The intensity is shown here instead of fields because the former is an experimental observable. Both the supermodes corresponding to the two eigenvalues show identical intensity, although the fields depict the bonding and antibonding character. See Supplementary information for the field distributions. With reducing gain, the splitting reduces proportionately. For example, the middle panel of (a) shows the situation for $n_{im} = -0.0073$ in cavity 2, where the $\delta\lambda = 1~$nm. The QE oscillation frequency and the phase difference ($\Delta T/ T$) diminishes. The intensity distribution for the two supermodes is skewed in favour of the cavity 1, as seen in (b).  Ultimately, at $n_{im} =  -0.0056$ in cavity 2, the system hits an exceptional point. The temporal emission coalesces into a single profile seen in the bottom panel of (a). The two eigenvalues coalesce and a single spectral peak is observed. The entire energy is accumulated in the cavity 1, which is indicative of eigenvector coalescence seen in the field distributions. See Supplementary Information for the fields. Fig\ref{fig:simring}(c) describes the systematic variation of the phase difference between two echoes and the output intensity with $\delta\lambda$.  A monotonic decrease in $\Delta T/T$ (blue markers, left Y-axis) with $\delta\lambda$ is seen. The peak output intensity $I_p$ increases upon approaching the EP, and maximises at the EP, indicative of enhanced functionality of the system operating at the exceptional point. The inset shows that $\Omega$ reduces monotonically with $\delta\lambda$. The correspondence of the quantum echo with the PT-symmetry breaking was confirmed by calculating the complex eigenvalues of the system. These  were computed by using a finite-element eigensolver computation, wherein the exact same structure was simulated. The gain/loss parameter $n_{im}$ was same as in the FDTD. The real (blue curve) and imaginary (red curve) eigenvalues are shown in subplot (d). The circle markers on the blue curve are the situations discussed in (a) and (b). A clear correspondence between the temporal behavior and the theoretically-expected PT-symmetry breaking and EP-formation is seen.

\begin{figure}[h!]
\centering
\includegraphics[width=0.9\textwidth]{1Dapprox_arxiv.pdf}
\caption{Schematic of linear arrays of spherical microcavities: (a) white arrows show emission and scattered directions: longitudinal emission (transmitted signal) and out-of-axis emission. The sizes of the cavities are about 18~$\mu m$. (b) A Finite-Element-Method (FEM) simulation showing the longitudinal mode.}
\label{fig:1Dapprox}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=1\textwidth]{Exp_setup_arxiv.pdf}
\caption{\textbf{Experimental setup and characterization of the Anderson localized lasing system:} (a) Complete schematic of the setup, see main text for details; (b) CCD image of the microresonator array; (c) zoomed view of the rectangular region in (a): the microresonator array is pumped by a mode-locked Nd:YAG laser (green beam, $\lambda = 532~nm$, pulse width = $30~ps$, repetition rate = $10$~Hz). The emission in the longitudinal direction (transmitted light) is used for the temporal measurement $I_{tmp}$, and the transverse out-of-axis scatter  gives the spatial profile $I_{sp}$ simultaneously. The schematic of emission directions is shown in (d)). The transverse scatter is imaged onto the input slit of a spectrometer S1 coupled with ICCD using a 4f imaging setup [L2: imaging lens  (10 cm focal length); F2: Notch filter (to remove the pump beam)]. The longitudinal emission is measured using a mirror M located at an angle of $45^{0}$ w.r.t the array as shown in the schematic. This emission is directed towards the input slit of a second identical spectrometer S2 coupled to a streak camera (Optronis SC-10, temporal resolution = $\sim3$~ps)  through a focusing lens L3 and a notch filter F2.  (e, f) represent the spatial data, where (e) shows the spectrographic image of the multiple lasing modes, and (f) shows a typical extracted spatial profile of one mode, yielding localization length; (g) shows the simultaneous streak image of the same modes with (h) showing a typical extracted temporal profile. (i) Lasing spectrum obtained by integrating the streak image over time axis, or integrating the spectrographic image over spatial axis. }
\label{fig:expsetup}
\end{figure}

The structures we employ for these experiments, schematized in Fig~\ref{fig:1Dapprox}(a), are linear arrays of microcavities that couple to form photonic bands and bandgaps. The microcavities are spherical in shape with diameters of about $18~\mu$m. The large diameters of the spheres and the linear configuration of the array induces coupling of longitudinal Fabry-Perot modes of the microcavities.  Figure~\ref{fig:1Dapprox}(b) illustrates a finite-element computation of an eigenmode of the periodic array, clearly depicting the longitudinal mode with quasi-planar wavefronts, akin to a multilayer system. Loss in this structure is induced by out-of-axis scattering of the mode by the interfaces. In order to introduce gain, we choose the  material appropriately to make the microcavities, as described in the Methods.  Disorder is introduced by deliberately modifying the axial dimensions and spatial separations of the microcavities. In our earlier works, we have shown that the periodic arrays sustain Bloch modes\cite{joshi2022anomalous}, while the disordered arrays realize Anderson localized modes\cite{joshi2020reduction, joshi2019effect, kumar2017temporal}. As the fabrication process of the microcavity arrays is dynamic in nature, we can monitor a huge number of configurations, as discussed next.

Figure~\ref{fig:expsetup} illustrates the creation of the microcavity array, and its simultaneous spectro-spatio-temporal characterization. The generation technique, as described in details in the Methods section, creates the microcavities from a Rhodamine-in-alcohol solution. The Rhodamine molecules act as sources distributed throughout the length of the array, whilst also providing gain to the mode. [CCD image shown in (b)]. Optical excitation is provided by pulses of an Nd:YAG laser fired at 1 Hz, illustrated by the green beam. The emission transverse to the array comprises the out-of-axis scatter that is imaged onto a identical spectrometer (S1) coupled to an imaging CCD. The longitudinal emission of the array, which carries the temporal signature of the emission in Anderson localized modes, is directed to an identical spectrometer (S2) coupled to a Streak Camera SC (Optronis SC-10).  Since the microcavities are constantly created and introduced into the pump laser focus, every excitation laser pulse sees a different disorder configuration of the array. This feature enables us to generate the vast ensemble of configurations required in this study. Importantly, the Streak Camera and Spectrometer CCD are synchronized with the pump laser to register the spectral, spatial and temporal data of every configuration simultaneously. Panels (c) and (d) illustrate the characterization scheme. Panel (e) depicts the spectrographic image from S1 at a particular pump pulse. It shows the multiple emission modes from the corresponding configuration illuminated by the laser pulse. The horizontal axis represents the spatial axis, while the  vertical axis represents the wavelength. Subplot (f) shows the extracted spatial profile of a representative mode labeled by the white arrow in the image (e). This is an Anderson-localized mode, whose localization length $\xi$ can be extracted from the decay in the wings. Panel (g) shows the Streak image of the emission of the same configuration. Since this is directed through a spectrometer, the  vertical axis  corresponds to the wavelength, while the  horizontal axis corresponds the time axis. A comparison between (g) and (e) shows that the same modes are captured in both devices, endorsing the simultaneous diagnostic capability of the setup. Subplot (h) is the extracted temporal profile of the same mode as in (f). Subplot (i) shows the emission spectrum obtained by integrating over the time-axis of (g). The same can be obtained by integrating the space-axis of (e).

We used this diagnostic system to characterize 3,000 different configurations of disorder, providing a massive ensemble of about 15,000 Anderson-localized lasing modes. The limit on the configuration number is set by the volume of the liquid present in the generator, as mentioned in the Methods. Although all modes were captured at the same pump energy, inherent pump fluctuations motivated us to further process the ensemble to choose modes in a very narrow range of pump energy, to within $\sim1~\mu$J at the average pump energy of $\sim87~\mu J$. The filtered ensemble still had over 5,000 modes. Among all the captured modes, pairs of coupled modes were identified through the temporal profiles, since the coupling resulted in a quantum echo. This clearly stood apart against the smooth profiles of individual uncoupled modes, such as the one shown in Figure \ref{fig:expsetup}(h). The coupling was further corroborated by the splitting in the spectral profile. Various degrees of splitting allowed us to identify the approach to a second-order exceptional point, and eventually the EP itself.

\begin{figure}[h!]
\centering
\includegraphics[width=1\textwidth]{experimental_sptemp_coupled_arxiv.pdf}
\caption{\textbf{Experimentally measured temporal, spectral, and spatial characteristics of Anderson localized lasing modes around and at an exceptional point (EP):} Subplots (a-d) show the measured quantum echo from lasing emission, with (d) being the EP. Unlike Figure ~3(h), the temporal profiles in (a) exhibit a quantum echo that decisively certifies that the two modes are coupled. Reducing echo frequency with decreasing phase difference between the echoes at the two wavelengths is clearly observed. Subplots (e-h) depict the corresponding eigenvalue spectra from the same configurations as (a-d) with mentioned splittings. The echo behavior is clearly consistent with the splitting. (i-l)  represent corresponding spatial intensity distribution of the coupled localized states. Each spatial mode exhibits two peaks, and signatures of shape exchange. Energy is distributed in both Anderson cavities at large splittings, but resides mostly in one cavity at $\delta\lambda = 0.2$~nm. A single profile exists at the EP in (l). The intensity behavior certifies coalescence of eigenvectors towards and at the EP. Subplots (d), (h) and (l) represent Anderson-localization lasing over an exceptional point, further confirmed in Figure 5.}
\label{fig:sptemp}
\end{figure}

Figure~\ref{fig:sptemp} illustrates a representative set of measurements in the Anderson-localizing system showing the approach to an EP. The first column of panels shows the measured quantum echo, the middle column shows the corresponding spectral measurement, while the third column illustrates the corresponding spatial intensity distributions. Each row of three panels show simultaneous diagnostics from the same disorder configuration, while different rows represent data from different configurations. The subplot (a) shows a temporal profile in which a high frequency quantum echo is observed. The two profiles (orange and blue) arising from two coupled modes are out of phase by $\Delta T/T$ as indicated in the plot. Subplot(e) shows the spectral splitting in the emission from the two coupled modes in this configuration. Quite symmetric spectral peaks are observed. The marked $\delta\lambda = 0.6$~nm was measured using a peakfinder algorithm. Subplot (i) shows the spatial intensity of the two modes extracted from the spectrographic data as described earlier. Strong spatial overlap with signatures of separated peaks are seen in the profiles. The coupled Anderson localized modes are spread over almost the system length. Two distinct peaks in both the profiles are seen, indicative of the two resonant cavities that underwent coupling. These profiles strongly emphasize the fact that the two coupled cavities are not physically separated, unlike conventional photonic coupled-cavities (Figure 1). In such a disordered system, the gain/loss values in the two cavities are self-determined by the random configuration that decides the relevant parameters such as quality factors, pump distributions etc. The lower rows show gradually decreasing quantum echo frequency, which corresponds to the reduction in splitting seen in the middle panels. At the same time, the intensity distributions morphologically coalesce as the splitting diminishes. Notably, subplot (k) shows a strong peak accompanied by a shoulder, that can be interpreted as the intensity preferentially residing in one cavity. These configurations represent the approach to a second-order exception point. Ultimately, the bottom row (d, h, l) represents lasing {\it at} the exceptional point, showing only one echoless temporal profile, one spectral peak and one intensity profile. Subplots (a) through (d) illustrate the vanishing of the echo as the EP is hit. The conclusive inference that (d) corresponds to a second-order EP is provided in the next figure where we present statistical analysis of the data.

\begin{figure}[h!]
\centering
\includegraphics[width=1\textwidth]{Lorentzfit_arxiv.pdf}
\caption{\textbf{Confirmation of the EP in a statistical ensemble:} (a) Scatter points for the QE frequency $\Omega$ show an overall reduction with decreasing $\delta\lambda$. Pink diamonds indicate $\langle \Omega \rangle$ at a particular $\delta\lambda$. (b) Measured $\Delta T/T$ and (c) peak emission intensity ($I_p$) as a function of $\delta\lambda$, measured over $\sim100$ modes. Each datapoint in (b) is an average over a few oscillations in the temporal profile. Pink squares indicate $\langle \Delta T/T \rangle$  at a particular $\delta\lambda$. The scatter clearly shows an approach to  $\Delta T/T \rightarrow 0$ as $\delta\lambda \rightarrow 0$. The yellow band around $\delta\lambda = 0$ signifies the EP region. All datapoints here are placed at $\delta\lambda = 0$ as the $\delta\lambda$ is beyond experimental spectral resolution. Although $\Delta T/T$ does not exist here, a suggestive datapoint is marked at $\Delta T/T = 0$ in logical follow-up to the neighbouring datapoints. (c) shows increasing peak intensity $I_p$ towards the EP region, with significant jump of over $300\%$ in the $I_{p}$ in the EP region. Inset shows $P(I_p)$ with emphasis shading on outlier points of high intensity. (d) Spectral line shape at the EP (filled green plot). Blue and red plots are Lorentzian and squared-Lorentzian fits respectively. The profile possesses a clear squared-Lorentzian character, which is the classic signature of an EP. Inset shows the same in logarithmic Y axis.}
\label{fig:lorentz}
\end{figure}

Fig~\ref{fig:lorentz}(a) illustrates the frequency of the quantum echo for various splittings. The mean frequency $\langle \Omega \rangle$ at each $\delta\lambda$ is marked as pink diamonds. The trend in $\langle \Omega \rangle$ is determined by the distribution of the scatter points which are constrained to a lower limit imposed by the total pulse duration ($\sim 50$~ps). (b) shows the phase difference between the echograms of the coupled modes as a function of $\delta\lambda$ over several configurations.  The scatter in the data reflects the probabilistic variation in the phase difference for a given splitting, with a large scatter existing in the range $0.2$~nm$<\delta\lambda<0.4$~nm, wherein a large number of coupled modes were observed. Such statistical variations are expected in non-deterministic systems due to configurational variations. The $\langle \Delta T/T \rangle$ is marked as pink squares. As $\delta\lambda \rightarrow 0$~nm, the $\langle \Delta T/T \rangle \rightarrow 0$. The $\langle\delta T/T\rangle$ shows a steady rise with $\delta\lambda$. The yellow band at $\delta\lambda = 0$ signifies the EP range. Clearly, $\delta T/T$ is undefined here since there is only one mode and one temporal profile. Nonetheless, a datapoint is intentionally marked to signify the completeness of the trend. Figure \ref{fig:lorentz}(c) shows the peak intensity of the modes. The intensity rises as modes begin to coalesce spectrally, in agreement with the typical behavior in an EP-approach. Inside the EP range (yellow band at $\delta\lambda = 0$), six modes were identified that lased at a very high intensity. Upto 300\% enhancement in lasing intensity was observed at these points. Such high intensity is an indicator of enhanced functionality, a preliminary identifier of a possible EP. The inset shows the intensity distribution over the entire ensemble of modes, emphasizing the outlier nature of the high-intensity points. Finally, subplot (d) confirms the exceptional point character of the modes using an analytical fit for the lineshape of one of the high-intensity peaks. It is known that LDOS is modified at the exceptional point such that the resonant lineshape approaches to the $n^{th}$-power of Lorentzian for an $n^{th}$-order EP. Consistent with this, the red solid line in (d) depicts the square-Lorentzian profile $\propto \frac{\gamma^2}{4\pi[(\Delta\omega^2 + \gamma^2)]^2}$ that perfectly fits the experimental peak. In comparison, the blue line is the conventional Lorentzian fit which visibly deviates from the observed data. The inset shows the same on a logarithmic Y-axis, emphasizing the fit by the square-Lorentzian. The modified LDOS at the EP leads to a maximal peak enhancement of 4 in passive systems, and can be even more in the presence of gain. Although our system is random and every configuration has cavities of varying quality-factors, the average behavior of $I_p$ in subplot (c) does exhibit instances of peak enhancement $>4$ compared to the coupled modes. The Supplementary material elucidates further data regarding fits on other exceptional point modes. Thus, these spectral lineshape diagnostics confirm that the lasing modes shown here occurred over exceptional points in non-Hermitian Anderson-localizing systems.

In summary, we have identified exceptional points in one-dimensional Anderson localizing systems, and have extracted lasing over the corresponding modes. The methodology adopted here significantly deviates from the tuning-parameter mechanism that is conventionally utilized in photonics. Our technique preserves the statistical variations in disorder that are sacrosanct in the domain of mesoscopic optics. We sample thousands of disorder configurations, and identify the quantum echo which is set up whenever coupled Anderson-localized modes are realized within the structures. The echo provides unambiguous signatures of coupling, and sets the coupled modes apart from random neighbouring modes in the spectrum. Our method realizes a large multitude of disorder configurations and simultaneously measures spatial, spectral and temporal dynamics therein. The presence of a quantum echo in the temporal signal allows us to track configurations with spectral splitting. We find that reduced splitting is accompanied by enhanced output intensity which allows us to identify individual modes that provide exceedingly large outcoupled intensity. Spectral shape analysis of these intense modes reveal the characteristic square-Lorentzian shapes originating from a second-order exceptional point degeneracy, as against Lorentzian peaks otherwise. The excess intensity in these modes is testimonial to the enhanced functionality afforded by the exceptional points.

\section*{Methods}

\subsection*{Sample preparation and system characterization}
The system comprises an array of amplifying micro-droplets made of Rhodamine-6G (with a concentration of 1.0~mM) dissolved in a mixture of methanol and ethylene glycol with an equal proportion(50:50$\%$). Ethylene glycol(EG) is added to obtain good stability and high index contrast (refractive index of the droplet, $n_{d}$ = 1.38) due to its high viscosity and refractive index value. Rh6G has a broadband gain profile in the range of 540$-$610~nm.

The array is generated using a vibrating orifice aerosol generator (VOAG). The liquid from a chamber (1000~ml liquid capacity) was pressurized using a dry-nitrogen cylinder and passed through VOAG. The VOAG has a circular opening with a diameter of $10~\mu m$, which generates an inhomogeneous and unstable cylindrical jet of droplets. In order to make stable and spherical-shaped droplets,
a piezoelectric gate is attached at the circular opening of the VOAG. The gate is then perturbed by applying a periodic electronic signal. Depending upon the amplitude and frequency of the signal, the size and spacing of the droplets can be varied in a controlled way which allows us to generate an ensemble of monodisperse droplets which acts as a periodic array. Controlled deviations in the electronic signal perturb the microdroplets and form a disordered array. For this work, we generate a disordered array of droplets. We made a large ensemble of disordered arrays with a mean diameter of 18.6~$\mu$m. The standard deviation in the diameters is 340.0~$nm$ which determines the strength of disorder in the array.

\subsection*{Coupled mode analysis of the binary system:}

1. Eigenvalue solver using a finite element method (FEM: COMSOL Multiphysics) was used to calculate the complex eigenfrequencies of the coupled ring resonator. Simulations were performed in 2D. The structural parameters of the coupled cavities are: Diameter of the ring $D = 2.0~\mu m$, thickness $t = 0.2~\mu m$ and separation $s = 0.4~\mu m$. The cavity resonance frequency for these parameters is $\sim$190.35~THz. The cavities' gain or loss are provided in terms
of their complex refractive indices, $n = n_{r}+ in_{im}$. The real part of the refractive index ($n_{r}$) for both cavities is set at 3.3. The separation between the cavities is also fixed at $s$, thus fixing the coupling strength ($\kappa$) between the cavities. We provide a fixed amount of gain ($n_{im} = -0.01$) in the $1^{st}$ cavity. The simulation was then performed varying $n_{im}$ of the second cavity from -0.01 to 0. During the tuning, exceptional points were identified to occur at $n_{im}\sim-0.0056$.

2. Temporal analyses were done using Finite-difference time-domain method (Lumerical FDTD). The structural and material parameters are same as mentioned in part 1 above. In this method, we use a Lorentz gain model with a Lorentz permittivity and resonance frequency close to the cavity frequency of the ring resonator. The permittivity of the material for this method is given by the relation:
\begin{equation}
\epsilon(\omega) = \epsilon + \epsilon_{l}\omega_{0}^{2}/(\omega_{0}^{2}-2i\delta_{0}\omega-\omega^{2})
\end{equation}

The gain center frequency is set at $\omega_{0}$ (Lorentz Resonance) while the width is set as $\delta_{0}$. The strength of the gain (gain amplitude) is given by the permittivity $\epsilon_{l}$, while the sign determines whether the material has loss or gain.

The gain level ($n_{l} = \sqrt{\epsilon_{l}}$) in the $1^{st}$ cavity is kept fixed at -0.01. It was confirmed that this value of gain did not lead to any divergence due to an exponential growth of the field. The gain level in the second cavity is gradually reduced from this value and replaced with loss. Both cavities are uniformly excited using a set of identical dipole sources at the cavity resonance frequency. The time signal of the field and intensity are measured by placing field time monitors inside the cavities. Spectra are calculated using the standard Fourier transform method. A full apodization method is used to measure the frequency-resolved temporal evolution of the cavity modes. This technique applies a window function to the fields $E(t)$ before its Fourier transform. By using this method, it is possible to calculate $E(\omega)$ from a portion of the time signal.

\bmhead{Data availability}
All data relating to the paper are available from the corresponding author upon reasonable request.

\bmhead{Code availability}
All the relevant computing codes used in this study are available from the
corresponding author upon reasonable request.

\bmhead{Supplementary information}
A Supplementary Information document has been attached with the manuscript.

\bmhead{Author contributions}
S.~M. conceived the project. S.~M. and K.~J. designed the experimental setup.  K.~J. performed the experiment, numerical simulations and analyzed the data.  K.~J. and S.~M. wrote the manuscript. S.~M. supervised the whole project.

\bmhead{Competing interests}
The authors declare no competing interests.

\bmhead{Acknowledgments}
SM would like to acknowledge the Swarnajayanti Fellowship from the Department of Science and Technology, Government of India. We acknowledge funding from the Department of Atomic Energy, Government of India (12-R $\&$ D-TFR-5.02-0200).