\section{Introduction}

One of the greatest successes in cosmology is the precise measurements of the cosmic microwave background (CMB) which support the inflationary-universe paradigm combined with the $\Lambda$-cold-dark-matter ($\Lambda$CDM) model. However, as data improves, a significant discrepancy emerges between its Hubble constant $H_0=67.36 \pm 0.54$ km/s/Mpc determined by Planck 2018 (P18)~\cite{Aghanim:2018eyx}, and $H_{0}=73.3 \pm 0.8$ km/s/Mpc obtained from the local distance-redshift (DR) measurements (see~\cite{Verde:2019ivm} and also~\cite{Knox:2019rjx,DiValentino:2021izs} and references therein). At the same time, the $\Lambda$CDM/P18 data fitting gives the clustering $S_8=0.832 \pm 0.013$~\cite{Aghanim:2018eyx}, while the recent weak-lensing (WL) data of KiDS-1000 and DES yield $S_8=0.766^{+0.020}_{-0.014}$~\cite{Heymans:2020gsg} and $0.773^{+0.026}_{-0.020}$~\cite{Abbott:2017wau}, respectively.

The axi-Higgs model recently proposed~\cite{Fung:2021wbz} suggests a potential solution to the $H_0$ tension by coupling ultralight axions to the Higgs field. This model can further explain the $^7$Li puzzle in Big-Bang nucleosynthesis (BBN), as well as the CMB $S_8$ tension with the WL data and the observed isotropic cosmic birefringence (ICB) in CMB~\cite{Minami:2020odp}. In this letter, we will demonstrate how this model resolves the $H_0$ and $S_8$ tensions simultaneously, by correlating the axion physics in the early and late universe.

For the application here, keeping all parameters in the standard model of strong and electroweak interactions unchanged, the axi-Higgs model with a single axion-like particle $a$ and the electroweak Higgs doublet $\phi$ is given by
\begin{eqnarray}\label{model1}
	V(a)  &=& m^2 a^2/2 + \left|m_s^2F(a) - \kappa \phi^{\dagger} \phi\right|^2 \ ,\\
F(a) &=& (1+\delta v)^2 = 1+ C{a^2}/{M_{\rm Pl}^2} \
\end{eqnarray}
where the axion mass $m$ is  $\sim 10^{-30} - 10^{-29}$~eV~\cite{Fung:2021wbz}; $m_s$ and $\kappa$ are fixed by the Higgs vacuum expectation value (VEV) today $v_0=246$ GeV and the Higgs boson mass $m_{\phi}=125$ GeV; $\delta v=\Delta v/v_0=(v-v_0)/v_0$ is the fractional shift of $v$ from its today's value $v_0$; $M_{\rm Pl}=2.4 \times 10^{18}$\,GeV is the reduced Planck mass.

Note that, because of the perfect square form of the Higgs potential, the electroweak contribution to the vacuum energy density is precisely zero, a necessary condition for an exponentially small cosmological constant without fine-tuning. This perfect square form is also crucial in the axi-Higgs model, where it suppresses the impact of the Higgs evolution on the axion, so $a$ evolves as if it is free. On the other hand, the presence of the axion (1) lifts $v$ in early universe and relaxes it to today's value in late universe, (2) contributes to the cold dark matter density which explains the ICB data, (3) because of its long wavelength, dampens the clustering amplitude, and (4) provides observable tests via atomic clock and/or quasar spectral measurements. Although it is known that $v > v_0$ may help to solve the $^7$Li puzzle in BBN \cite{Kneller:2003xf, Li:2005km,Coc:2006sx,Dent:2007zu, Browder:2008em,Bedaque:2010hr,Cheoun:2011yn,Berengut:2013nh,Hall:2014dfa,Heffernan:2017hwa,Mori:2019cfo} and the Hubble tension \cite{Ade:2014zfo,Hart:2019dxi}, the axi-Higgs model provides the mechanism how $v>v_0$ happens before recombination and how $v$ relaxes to $v_0$ today.

The axion field $a$ stays in a misaligned state $a_{\text{ini}}$ initially until the Hubble parameter drops to $H(z) \sim m$. Then $a$ rolls down along the potential, and deposits its vacuum energy into dark matter (DM). This happens at the redshift $z \sim 1000 - 100$, and determines today's relic abundance $\omega_a$. Since $z_\text{BBN} \sim 10^9$ at BBN while $z_* \sim 1100$ at recombination, $a_{\text{ini}}$ contributes to $\delta v_{\text{ini}}=\delta v_{\text{BBN}}=\delta v_{\text{rec}}$ from the BBN epoch to the recombination epoch, where $\delta v_\text{BBN} \sim 0.01$ can largely solve the $^7$Li puzzle~\cite{Fung:2021wbz}. Replacing the coefficient $C$ and $a_{\text{ini}}$ by $\delta v_{\text{rec}}$ and $\omega_a$ ($C \simeq 0.12 \delta v_{\text{ini}}/\omega_a$), we have five relevant parameters to be determined: the baryon density $\omega_b$, the cold dark matter density $\omega_c$, $h=H_0/100$~km/s/Mpc, $\delta v_\text{rec} (=\delta v \ {\rm below})$ and $\omega_a$, with
\begin{equation}
h^2 = \sum \omega_i = \omega_b + \omega_c + \omega_a + \omega_\Lambda \ .   \label{eq:freq}
\end{equation}

    \begin{table}[htp]
\resizebox{0.47\textwidth}{!}{
			\begin{tabular}{|c|cc|cc|}
				\hline
				\multirow{2}{*}{} & $\Lambda$CDM & axi-Higgs & $\Lambda$CDM & axi-Higgs \\
				& \multicolumn{2}{c|}{(CMB+BAO+WL+SN)} & \multicolumn{2}{c|}{(CMB+BAO+WL+DR)} \\
				\hline\hline
				$\omega_b$ & $0.02252 \pm 0.00013$ & $0.02276 \pm 0.00021$ & $0.02267 \pm 0.00013$ & $0.02299 \pm 0.00014$ \\

				$\omega_c$ & $0.11802 \pm 0.00080$ & $0.1209 \pm 0.0019$ & $0.11657 \pm 0.00077$ & $0.1228 \pm 0.0014$ \\

				$H_0$ & $68.24 \pm 0.36$ & $70.3 \pm 1.6$ & $68.96 \pm 0.34$ & $72.42 \pm 0.76$ \\

				$v / v_0$ & $1$ & $1.0142 \pm 0.0088$ & $1$ & $1.0254 \pm 0.0050$ \\

				$1000\omega_a$ & $0$ & $< 1.28~[95\%]$ & $0$ & $< 1.32~[95\%]$ \\
				\hline

				$\omega_\Lambda$ & $0.3244 \pm 0.0056$ & $0.349 \pm 0.021$ & $0.3356 \pm 0.0053$ & $0.377 \pm 0.010$ \\

				$\omega_m$ & $0.14117 \pm 0.00077$ & $0.1448 \pm 0.0022$ & $0.13988 \pm 0.00075$ & $0.1470 \pm 0.0017$ \\

				$S_8$ & $0.8082 \pm 0.0092
				$ & $0.8035 \pm 0.0097$ & $0.7902 \pm 0.0088$ & $0.7970 \pm 0.0088$ \\

				$100\theta_*$ & $1.04129 \pm 0.00029$ & $1.04114 \pm 0.00030$ & $1.04154 \pm 0.00028$ & $1.04107 \pm 0.00030$ \\

				\hline
				$-\ln \mathcal{L}$ & $8.4$ & $5.9$ & $30.4$ & $16.2$ \\
				$N_Y - N_X$ & $22 - 3$ & $22 - 5$ & $23 - 3$ & $23 - 5$ \\
				\hline

			\end{tabular}
}
		\caption{Marginalized parameter values in the $\Lambda$CDM model and the axi-Higgs model with $N_Y$ is the number of observables included and $N_X$ is the number of parameters involved. The upper bounds of $\omega_a$ are shown with 95\% confidence level.}
		\label{tab:model_comparison}
	\end{table}

Considering that the parameters determined in the axi-Higgs model deviate only slightly (at a percent level) from those in the $\Lambda$CDM model, we apply a leading-order perturbative approach (LPA)~\cite{Fung:2021wbz} to analyze the Hubble flow, i.e., the evolution of $H(z)$ as a function of $z$. This semi-analytic approach allows one to see clearly the inter-dependence of the model parameters versus the observational data. A full data analysis by one of us ~\cite{Luu:2021} provides a valuable check on the validity of the LPA.

As to be seen, the impacts of varying the other parameters beyond Eq. (\ref{eq:freq}) are sub-leading. So we simply fix them to the default/best-fit values in $\Lambda$CDM/P18~\cite{Aghanim:2018eyx}, which include the density $\omega_{\nu}$ for two massless neutrinos and one light one with $m_{\nu} = 0.06\,$eV; $10^9A_s = 2.10055$ and $n_s = 0.96605$ to the initial curvature spectrum, and the reionization optical depth $\tau_\text{re} = 0.05431$. The flatness of the universe is assumed throughout the universe history.

The analysis results are summarized in Tab.~\ref{tab:model_comparison} and Fig \ref{fig:model_comparison} with $m=2 \times 10^{-30}$ eV as a benchmark, where the axi-Higgs model has considerably improved likelihoods compared to the $\Lambda$CDM model. Combining the CMB+BAO (baryon acoustic oscillation)+WL+SN(supernovae) data yields the axi-Higgs best fit
$H_0 = 70.3 \pm 1.6$~km/s/Mpc and $S_8 = 0.8035 \pm 0.0097$. Especially, here $\omega_a \simeq 0.0004$ is sufficient to explain the ICB data,

 and $\delta v=1.42$\% is in good agreement with $\delta v_{\text{BBN}}=1.2$\% required to solve the $^7$Li puzzle \cite{Fung:2021wbz}. However, the DR data further up-shifts $H_0$ with $\delta v = 2.54$\%, which  is higher than needed by BBN.
 If this $\delta v$ tension persists, it can be solved with the introduction of another axion~\cite{Fung:2021wbz}, conveniently the one ($m \simeq 10^{-22}$ eV) for fuzzy dark matter~\cite{Hu:2000ke,Schive:2014dra,Hui:2016ltb}.

	\section{The LPA Analysis}

With the set of $N_X$ parameters characterizing the model and a collection of $N_Y$ ($\ge N_X$) observables representing the data, we derive variation equations of these observables in terms of the set of parameters, with respect to (w.r.t.) a specific reference point. We choose the best fit in the $\Lambda$CDM/P18 scenario~\cite{Aghanim:2018eyx} as the reference point. The equations are then given by
		\begin{eqnarray}
			\delta Y = Y_{|b} \delta \omega_b + Y_{|c} \delta \omega_c + Y_{|h} \delta h + Y_{|v} \delta v + Y_{|a} \omega_a \; .  \label{eq:gen}
		\end{eqnarray}
		Here $\delta Y \equiv d\ln Y$ and $Y_{|X} \equiv \partial\ln Y / \partial\ln X$, with the exception of $Y_{|a} \equiv \partial\ln Y / \partial \omega_a$. These coefficients  are calculated either analytically from their definitions or numerically using public Boltzmann codes. The $Y_{|X}$ values used in this study are summarized in Tab.~\ref{tab:par_derivs_final} (see supplemental material at Sec.~A for a full list of $X$ and at Sec.~B for a full table of $Y_{|X}$).

	\begin{table}[htp]
\resizebox{0.45\textwidth}{!}{
		\centering
		\begin{tabular}{ |c|c|c|c|c|c| }
		\hline
		\diagbox{Y}{X}  & $ \omega_b$ & $ \omega_c$ &  $ h$ & $ v$ & $\omega_a$  \\
		\hline
	    $ l_a$ & $0.0550$ & $- 0.1203$ & $- 0.1934$ &  $0.6837$ & $- 2.6364$ \\
		\hline
		$ l_\text{eq}$ & $0.0942$ & $0.5082$ & $-0.1934$ & $0.0154$ & $- 2.6376$   \\
		\hline
		$ l_D$ & $0.2459$ & $- 0.0962$ & $- 0.1934$ & $0.4553$ & $-2.6315$ \\
		\hline
		$S_L$ & $-0.1398$ & $0.8484$ & $-0.2619$ & $0.0788$  &  $-16.082$  \\
		\hline
		$ \alpha_\perp$(0.698) & $0.1466$ & $0.0969$ & $-0.7175$ &  \multirow{3}{*}{$0.6162$} & $-0.9852$  \\
		\cline{1-4}\cline{6-6}
		$ \alpha_\parallel$(0.698) & $0.1256$ &$-0.0160$ & $-0.4483$ &  & $-1.9252$  \\
		\cline{1-4}\cline{6-6}
		$ \alpha_V$(0.845) & $0.1354$ &  $0.0370$ & $-0.5747$ &  & $-1.4838$  \\
		\hline
		$ S_8$ & $-0.1007$ & $1.0578$ & $-0.7658$ & $0.0788$  &  $-14.337$  \\
		\hline
		$ \mu (1.382)$ & $-0.0017$ & $-0.0089$ & $-0.0270$ & $0$  &  $-0.0743$  \\
		\hline
		\end{tabular}
		}
		\caption{$Y_{|X}$ values in the axi-Higgs model. The $\alpha(z_{\rm eff})$'s are for BAO and $\mu(z_{\rm eff})$ is distance modulus for supernovae. The $\Lambda$CDM(+$m_e$) model shares the values of  $Y_{|b,c,h}(+Y_{|v})$. }
		\label{tab:par_derivs_final}
	\end{table}

Then we apply a likelihood method to the group of variation equations, to find out the parameter values favored by data. The likelihood function is defined as
	\begin{equation}
		\mathcal{L} =
		\frac{1}{\sqrt{(2\pi)^{N_Y} | \mathbf{\Sigma}|}} \exp  \left[ -\dfrac{1}{2} ({\bf Y}_\text{o} - {\bf Y}_{\rm t})^\text{T} \mathbf{\Sigma}^{-1} ({\bf Y}_\text{o} - {\bf Y}_{\rm t}) \right]  \ .
		\label{eq:lik}
	\end{equation}
Here ${\bf Y}$ denotes the $N_Y$ observables, and its subscripts ``o'' and ``t'' represent observation values and model predictions, respectively. $\mathbf{\Sigma}$ is covariance matrix, given by $\mathbf{\Sigma}_{ij} = \rho_{ij} \sigma_i\sigma_j $, where $\rho_{ij}$ measures the observable correlation with $\rho_{ij} = \rho_{ji}$ and $\rho_{ii} = 1$, and $\sigma_i$ is observation variance. The numerical MCMC sampler \textbf{Cobaya}~\cite{Torrado:2020dgo} is used to sample the likelihood function for the analysis.

	 For the CMB data, we consider the sound horizon at recombination $l_a$, the Hubble horizon at matter-radiation equality $l_{\rm eq}$, the damping scale at recombination $l_D$ and $S_L$. $l_a$ determines the position of the first acoustic peak and also the peak-spacings. $l_\text{eq}$ sets up the threshold for radiation to dominantly drive the gravitational potential, while $l_D$ is the scale below which fluctuations are suppressed by photon-baryon coupling and multipole anisotropic stress. Phenomenologically, $l_\text{eq}$ and $l_D$ determine the relative peak heights while $l_\text{eq}$ also determines the modulation between the even and odd peaks. As pointed out in~\cite{Hu:1995kot, Hu:2000ti, Hu:2001bc}, the CMB power spectrum is essentially characterized by these three scale parameters.

Furthermore, in order to address the CMB constraints on matter fluctuation, we need to consider observables beyond the $C^{\rm{TT}}_{\ell}$ spectrum,
which is essentially encapsulated in the CMB lensing power spectrum $C^{\rm{\phi\phi}}_{\ell}$ and manifested as a combination $S_L \equiv \sigma_8 \Omega_m^{0.25}$~\cite{Planck:2018lbu}.

		\begin{table}[htp]
\resizebox{0.48\textwidth}{!}{
		\centering
		\begin{tabular}{ |c|c|c|c|c|c|c| }
			\hline
			${\bf Y}$ & ${\bf Y}_{\rm ref}$ & ${\bf Y}_{\rm o}$ & $\mathbf{\Sigma}$ & ${\bf X}$ & ${\bf X}_{\rm P18}$ & ${\bf X}_{\rm LPA}$\\
			\hline
			$l _a$ & $301.765$ & $301.757$ & $\multirow{4}{*}{Eq.~\eqref{Eq:cmb_cov}}$ & $\omega_b$ &$0.02237\pm 0.00015$ & $0.02237\pm 0.00014$ \\
			\cline{1-3} \cline{5-7}
			 $l_\text{eq}$ & $144.14$ & $144.05$ & & $\omega_c$ &$0.1200\pm 0.0012$ & $0.1200\pm 0.0011$ \\
			\cline{1-3} \cline{5-7}
			$l_D$ & $1954.6$ & $1954.3$ & & $h$ &$0.6736\pm 0.0054$ & $0.6735 \pm 0.0050$  \\
			\cline{1-3} \cline{5-7}
			$S_L$ & $0.6085$ & $0.6078$ & & $S_8$ & $0.832 \pm 0.013$ & $0.832 \pm 0.013$ \\
			\hline
		\end{tabular}
		}
		\caption{Test of the LPA validity in the $\Lambda$CDM model, with the best fit values of P18.}
		\label{tab:CMBdata}
		\end{table}

The variations of ${\bf Y} = \{l_a, l_\text{eq}, l_D, S_L \}$ can be read out directly, using the $Y_{|X}$ values in Tab.~\ref{tab:par_derivs_final} as the inputs of Eq.~(\ref{eq:gen}). Thus we have
\begin{eqnarray}\label{eq:dleq}
 \delta l_a  &=& l_{a|b} \delta\omega_b + l_{a|c} \delta\omega_c + l_{a|h} \delta h + l_{a|v} \delta v + l_{a|a}\omega_a \\
	&\simeq& 0.05 \delta\omega_b - 0.12 \delta\omega_c - 0.19\delta h +0.68 \delta v - 2.6 \omega_a \nonumber\\
 \delta l_{\rm eq} &\simeq& 0.09 \delta\omega_b + 0.51 \delta\omega_c - 0.19 \delta h + 0.02 \delta v - 2.6 \omega_a \nonumber\\
  \delta l_D &\simeq& 0.25 \delta\omega_b - 0.1 \delta\omega_c - 0.19 \delta h + 0.46 \delta v - 2.6 \omega_a \nonumber\\
\delta S_L &\simeq& -0.14 \delta \omega_b + 0.85 \delta \omega_c - 0.26 \delta h + 0.08 \delta v - 16 \omega_a \nonumber
\end{eqnarray}
and similarly for the other observables.
  Since $l_a$ is well measured with a very small uncertainty, a shift in $h$ has to be compensated by corresponding shifts in the other quantities to keep $\delta l_a\simeq 0$. On the other hand, with a large negative coefficient for $\omega_a$ in $\delta S_L$, lowering $\sigma_8$ needs only a small $\omega_a>0$.  The model predictions are then given by ${\bf Y}_{\rm t}= {\bf Y}_{\text{ref}} \left( 1 + \delta {\bf Y} \right)$. $\mathbf{\Sigma}$ can be calculated using the inputs from Tab.~\ref{tab:CMBdata} (see supplemental material at Sec.~C for details).

 We first validate our LPA approach by comparing the $\Lambda$CDM fitting results against P18's full analysis as shown in Fig.~\ref{fig:LPA}.
 The axi-Higgs model reduces to the $\Lambda$CDM model by simply setting $\delta v = \omega_a = 0$.
One can clearly see from this figure that the LPA exceptionally reproduces the marginalized contours of $\omega_b$, $\omega_c$ and $h$ and their posterior distributions from the P18's original fit~\cite{Aghanim:2018eyx}. The numerical LPA results (${\bf X}_{\rm LPA}$) are only slightly different from the best-fit of $\Lambda$CDM to P18 (${\bf X}_{\rm P18}$), in terms of both central values and their uncertainties (see Tab.~\ref{tab:CMBdata}). This justifies the LPA validity with only four observational inputs. We also note that dropping the correlations in $\mathbf{\Sigma}$ only causes a small change to the LPA contours in Fig.~\ref{fig:LPA}. This indicates that $l_a$, $l_\text{eq}$, $l_D$, $S_L$ as the CMB observables are relatively independent in $\Lambda$CDM. The LPA is equally successful when extended to the $\Lambda$CDM$+m_e$ model~\cite{Ade:2014zfo,Hart:2019dxi}, where $\delta v=\delta m_e$ (see supplemental material at Sec.~C for details). Since this method can be applied to simple extensions of the $\Lambda$CDM model with significantly less numerical effort, we will apply it to study the axi-Higgs model below.

	\begin{figure} [!ht]
		\centering
		\includegraphics[scale=0.48]{fig/LCDM_CMB.pdf}
		\caption{The LPA contours overlap with the P18/$\Lambda$CDM best-fit results~\cite{Aghanim:2018eyx} precisely. All posterior distributions of the model parameters in this letter are plotted with \textbf{GetDist}~\cite{Lewis:2019xzd}.}
		\label{fig:LPA}
	\end{figure}

\section{$\mathbf{H_0}$ and $\mathbf{S_8}$ in the Axi-Higgs Universe}

Despite the sharing of $m_e$ as a free parameter in the recombination epoch, the axi-Higgs model is essentially different from the $\Lambda$CDM+$m_e$ model, due to the impacts of the time-varying state of $a$, that is, as dark energy (DE) in the early universe and as DM in the late universe. For example, the $a$ evolution causes a positive shift to $v$ during $z > z_a$ and then brings it back to its today's value later~\cite{Fung:2021wbz}. In contrast, such a mechanism is lacking for the $\Lambda$CDM+$m_e$ model~\cite{Ade:2014zfo,Hart:2019dxi}. These impacts raise the $H_0$ value to a level higher than what the $\Lambda$CDM+$m_e$ model allows. The Hubble flow here is determined as
\begin{eqnarray}
H(z) &=& 100 \text{ km/s/Mpc} \times \\
 &&\sqrt{\omega_r (1+z)^4 + (\omega_c+\omega_b) (1+z)^3 + g(z) \omega_a + \omega_{\Lambda}} \ ,  \nonumber
\end{eqnarray}
where $\omega_r$ is the radiation density, and
\begin{equation}\label{eq:a1osc}
  g(z) = \begin{dcases}
    (1+z)^3 \ ,  & z \ll z_a<z_*  \ ; \\
   (1+z_a)^3 \ ,  & \quad  z \ge z_a  \ .
  \end{dcases}
\end{equation}

	\begin{figure}[!ht]
		\centering
		\includegraphics[scale=0.55]{fig/derivatives_plot.pdf}
		\caption{Evolution of $H(z)_{|a,b,c,\Lambda}$ and $\delta H(z)$ in the axi-Higgs benchmark, with CMB+BAO+WL+SN data.}
		\label{fig:phpw}
	\end{figure}
To reveal this mechanism, we display the evolution of $\delta H(z)=d \ln H(z)$ and its derivatives w.r.t. $\omega_{b, c, a,  \Lambda}$ in the axi-Higgs benchmark in Fig.~\ref{fig:phpw}, in which the axion oscillating feature that is ignored in the above Eq.(\ref{eq:a1osc}) is restored. One can see that the Hubble flow deviates from its $\Lambda$CDM prediction since $z>z_{\rm eq}$, which is sequentially taken over by $\omega_c$, $\omega_a$ and $\omega_\Lambda$ until today. This deviation is originally driven by $\omega_a$,
where $H(z)_{|a}$ evolves approximately for three stages. In the early time, the axion is DE-like. $H(z)_{|a}$ evolves as
$H(z)_{|a} \propto 1/ \omega_r (1+z)^4$ before the matter-radiation equality time and $\propto 1/(\omega_b +\omega_c)(1+z)^3$ after that, so its value is suppressed at high redshift. This lasts a long period until the axion becomes DM-like at $z \sim z_a$. $H(z)_{|a}$ then evolves roughly as a constant $\propto 1/\omega_m$ during $z \sim 100 - 1$, with a wiggling feature being developed for its curve due to axion oscillation. In the $\Lambda$-dominant epoch ($z < 1$), $H(z)_{|a}$ drops quickly as $z$ goes to zero, as $H(z)_{|a} \propto {(1+z)^3}/ {\omega_\Lambda}$. Such an evolution pattern of $H(z)_{|a}$, particularly its relatively big value after recombination, results in a universal negative dependence of the CMB and BAO scale parameters on $\omega_a$ (see Tab.~\ref{tab:par_derivs_final}).

Let us consider $l_{\rm eq} \propto H(z_{\rm eq}) D_*$ as an example. Here $H(z_{\rm eq})$ is determined by the early-time cosmology and hence less influenced by $\omega_a$, while the diameter distance $D_* =\int_{0}^{z_*} \frac{dz'}{H(z')}$ is closely related to the cosmic evolution after recombination, where the variation of $D_*$ caused by $\omega_a$ is $\sim \int_{0}^{z_*} \frac{ - H(z')_{|a} \omega_a dz'}{H(z')}$. This necessarily yields $l_{{\rm eq}|a} < 0$. Similar arguments can be applied to the other cases.

In the $\Lambda$CDM model, a value of $h$ or $H_0$ value close to its local DR measurements is highly disfavored by the CMB data because of its correlation with $\omega_b$ and $\omega_c$. The Friedman equation for the universe today (see Eq.~(\ref{eq:freq})) indicates that, as $h$ increases, $\omega_i$ tends to increase even faster. Being out-of-phase between these two variations breaks the variation equations of the CMB and BAO scale parameters defined by Tab.~\ref{tab:par_derivs_final}. This also worsens the existing $S_8$ tension due to its high sensitivity to $\omega_c$. However, the situation is different in the axi-Higgs model.

Among the CMB and BAO scale parameters, the said effect of varying $h$ yields the largest $|Y_{|c}|$ value for $l_{\rm eq}$. In the axi-Higgs model, this impact is significantly cancelled by the $\omega_a$ contribution. As expected from the discussion above, $\delta l_{\rm eq}$ (\ref{eq:dleq}) has an opposite dependence on $\delta \omega_{c}$ and $\omega_a$.
As for the impacts brought in by the requested $\omega_a$ on the other scale parameters, they will be largely absorbed by a positive $\delta v$ which also compensates for the impacts of shifting $h$ considerably. Except $l_{\rm eq}$, these parameters demonstrate a positive and comparable dependence on $\delta v$. Such a dependence mainly arises from the impacts of $\delta v$ on the sound horizon at recombination and the end of baryon drag, while $\delta v$ itself is driven by the $a$ evolution before $z>z_a$. By correlating axion physics in the early and late universe, the $H_0$ tension is thus reduced.

		\begin{figure} [!ht]
		\centering
		\includegraphics[scale=0.285]{fig/AxiHiggs_CMB-BAO-WL-SN.pdf}
		\caption{Posterior distributions of the model parameters in the axi-Higgs benchmark. The shaded olive and grey bands represent the local DR measurement of $H_0 = 73.2 \pm 1.3$ km/s/Mpc from the latest SH0ES-20~\cite{Riess:2020fzl} and the weak-lensing measurement of $S_8 = 0.766^{+0.020}_{-0.014}$ from KiDS-1000~\cite{Heymans:2020gsg}, respectively.}
		\label{fig:model_comparison}
	\end{figure}

	\begin{figure}[!ht]
		\centering
		\includegraphics[scale=0.4]{fig/axion_mass.pdf}
		\caption{Relative variations of the model parameters as a function of the axion mass $m$, in the axi-Higgs/CMB+BAO+WL +SN scenario. The central boxes and error bars denote the uncertainties at 68\% and 95\% significance level (see supplemental material at Sec.~D for more information).}
		\label{fig:axion_mass}
	\end{figure}

We demonstrate the posterior distributions of the model parameters in the axi-Higgs benchmark in Fig.~\ref{fig:model_comparison}. As a reference, we also show in this figure the LPA results of the $\Lambda$CDM and $\Lambda$CDM$+m_e$ models, with the given data.  Compared to the $\Lambda$CDM/CMB result, the $H_0$ value favored by the $\Lambda$CDM$+m_e$/CMB$+$BAO model is higher. This is consistent with~\cite{Ade:2014zfo,Hart:2019dxi}. However, the axi-Higgs/CMB$+$BAO$+$WL+SN model yields an up-shift of $H_0$ which is bigger than that achieved in the $\Lambda$CDM$+m_e$/CMB$+$BAO model, accompanied by a downward shift of $S_8$ (due to $S_{8|a} \ll 0$ and $\omega_a > 0$). Thus the blue filled contours overlap with the intersection of the shaded olive and grey bands in the $H_0$-$S_8$ panel - the $H_0$ and $S_8$ tensions are simultaneously reduced.

Actually, P18~\cite{Aghanim:2018eyx} measures $S_L$ while their extracted $S_8=S_L(\omega_m/0.30)^{1/4}/h^{1/2}$ value
uses the $h$ value determined by the $\Lambda$CDM model. The higher $h$ value determined in the axi-Higgs model lowers the $S_8$ value in CMB, thus weakening the $S_8$ tension before any other impacts from the axi-Higgs model are considered. But, the uncertainty introduced by the ``consistency parameter'' $A_L$ \cite{Aghanim:2018eyx} suggests that the full resolution of the rather weak $S_8$ tension may have to wait for a better understanding of clustering and weak lensing physics.

At last, we consider the impacts of the axion mass $m$, by varying its value from $10^{-28}$~eV to $10^{-31}$~eV. For $m = 10^{-28}$\,eV, we have $z_a\sim 2100$ and $z_{\rm eq} > z_a > z_*$. The impact of $a$ on $l_{\rm eq}$ is still given by Eq.~(\ref{eq:dleq}) approximately, but it starts to roll down before $z_*$. So keeping $\delta v_{\rm rec}$ at a few percent level requires a larger $C$ value. For $m =10^{-27}\,$eV, $a$ contributes to DM at $z_*$. At the smaller $m$ end, for $m =10^{-31}\,$eV, $a$ behaves like DE as its dropping is late, which results in a weakened suppression for galactic clustering or $S_8$. Notably, such a small $m$ has been ruled out by the AC measurements~\cite{Fung:2021wbz}.

\section{Summary and Remarks}

As a low-energy effective theory motivated by string theory, the axi-Higgs model broadly impact our understanding of the universe~\cite{Fung:2021wbz}. In this letter, we demonstrate how the $H_0$ and $S_8$ tensions get simultaneously relaxed in this model, by correlating the impacts of the axion $a$ on the early and late universe:

In the early universe ($z>z_a$), this axion field behaves like DE. Its main impact is to drive a positive shift of the Higgs VEV (see Eq.~(\ref{model1})). Such a shift will reduce the cross section of Thomson scattering ($\propto m_e^{-2}$) and modify various atomic processes crucial to electron recombination and photoionization as discussed in~\cite{Ade:2014zfo, Hart:2017ndk}, and therefore, e.g., increase $z_*$ and decrease the comoving sound horizon and damping scale at $z_*$.

 In the late universe ($z<z_a$),  this axion field behaves like DM. Its main impacts are to: (1) increase the $H(z)$ value during this epoch by providing additional contribution to the energy density and hence, e.g., reduce the comoving diameter distance at $z_*$ and $z_{\rm eff}$; (2) suppress the formation of the structure at a galactic clustering scale and even above, with its super-long de Broglie wavelength; (3) shift the $v$ (or $m_e$) value in the early universe to its today's value $v_0$; (4) generate the observed ICB~\cite{Minami:2020odp}.

These features of the axi-Higgs model distinguish itself from the class of popular proposals to address the $H_0$ tension, namely increasing the $H(z)$ values in the early universe with additional energy density contributions, such as early DE~\cite{Karwal:2016vyq,Poulin:2018cxd} and self-interacting neutrinos~\cite{Kreisch:2019yzn}.
To conclude, we emphasize that a full test of this model is at hand, due to the oncoming AC and the quasar spectral measurements with the data expected to be collected by, e.g., Thirty Meter Telescope~\cite{Skidmore:2015lga} and James Webb Space Telescope~\cite{Behroozi:2020jhj}. More details on this can be found in~\cite{Fung:2021wbz}.

\section*{Acknowledgement}

	We thank Luke Hart and Jens Chluba for valuable communications.
	This work is supported partly by the Area of Excellence under the Grant No.~AoE/P-404/18-3(6), partly by the General Research Fund under Grant No.~16305219, and partly by the Collaborative Research Fund under the Grant No.~C7015-19G.  All grants were issued by the Research Grants Council of Hong Kong S.A.R.

\section{SUPPLEMENTAL MATERIALS}

  In the brane-world scenario in string theory, this axion $a$ comes from the closed string sector, where the graviton also comes from. They are expected to be components of the complex-structure moduli responsible for the shape of the Calabi-Yau orientifold (and/or the dilaton mode). This is different from an axion in the open string sector inside a stack of (anti-D3-)branes, in which all standard model particles reside. In this perspective, we expect that the decay constant $f_a$ to lie somewhere between the GUT-string scale and $M_{\rm Pl}$. With $\omega_a/\omega_m \simeq 0.3\%$ (or $\omega_a \sim 0.0004$), addressing the ICB anomaly requests $a_{\text{ini}}/f_a \sim 1$ or $f_a \sim 2 \times 10^{17}$ GeV which is reasonable.


\section{Introduction}
	It has now been 60 years since John Stewart Bell wrote his famous paper on the {Einstein}-Podolsky-Rosen (EPR) paradox \cite{bell1964einstein}, and 50 years since the first experimental Bell test \cite{freedman1972experimental}. The majority of physicists are perfectly happy to concede that in the lab we see experimental results consistent with the postulates of quantum mechanics. However, the implications of these mathematical postulates on the `reality' of the wavefunction is still very much up for debate \cite{caves2002quantum,harrigan2010einstein,pusey2012reality,colbeck2012system,mermin2014physics,ringbauer2015measurements}.

	These Bell experiments remain as some of the most important demonstrations for the reality of the quantum state and the death of a `local realism' picture of nature. In such an experiment a physical system is distributed between spatially separated observers, and we allow these observers to perform measurements on their local system. The emerging statistics prove that physical systems are not bound to behave locally (in accordance to local hidden-variable models). Rather, the statistics are consistent with the postulates governing quantum {mechanics}.

	In this work we remove the postulates of quantum mechanics and instead allow a physical system to be distributed according to a {quasiprobability (hidden-variable) distribution that is allowed to take negative values}. Although we are perfectly content with real negative numbers in physics, negative {quasi}probabilities-- despite receiving support from individuals such as Dirac \cite{dirac1942bakerian} and Feynman \cite{feynman1987negative} and having a solid mathematical foundation \cite{ruzsa1988algebraic, khrennikov2007generalized}-- have been a long debated issue in theoretical physics \cite{muckenheim1986review}. See, for example, the extensive discussion surrounding the interpretation of negative values in the Wigner distribution \cite{ferrie2011quasi,veitch2012negative}. {{ In the majority of considerations, quasiprobability distributions are used to describe states that are not directly observed}; that is, all observable measurement statistics must be governed by ordinary probability distributions. {As an example{,} a Wigner function may assign a negative quasiprobability to a particle having a particular position/momentum combination, but any physical measurement, constrained by Heisenberg uncertainty, will have an all-positive outcome distribution.} This feature ensures that no outcome is ever predicted to be seen occurring a negative number of times \cite{feynman1987negative}, and similarly protects the quasiprobability physicist from falling victim to `Dutch book' arguments \cite[Ch.3]{de2017theory}}.

	An important motivator for this work is the result of Al-Safi and Short \cite{al2013simulating} (expanded upon by the authors of \cite{oas2014exploring}) which showed that it is possible to simulate all non-signalling correlations, (those which adhere to the principles of {special }relativity) \cite{popescu1994quantum,peres2004quantum}, if one allows negative values in a probability distribution. However, physical reality does not explore this full set of correlations -- but rather, is restricted to those achievable by quantum correlations. Therefore the question that we pose in this work is:\begin{quote}\center
		``\textit{What are the restrictions on the negativity in a {hidden-variable probability distribution} such that it can emulate the statistics seen in a physical Bell experiment?}''
	\end{quote}  \vspace{1em}
	In order to answer this question we  {construct
	CHSH inequalities} for two parties \cite{clauser1969proposed} whose degree of violation is witnessed by the amount of negativity present in the hidden-variable probability distribution. Our witness { yields}  a value of 0 for a quasiprobability distribution which is entirely positive, such as that which
	would describe an ordinary classical system.

	We start by describing the setup necessary for the construction of these nonlocal experiments, introducing the probability distributions admitted by classical, quantum, and non-signalling {theories}, and giving the famous Bell scores that these distributions respectively allow one to reach in such nonlocal experiments.  We then give the definition of a quasiprobability distribution and motivate \textit{negativity witnesses} as a quantitative method of detecting negativity in said distributions. Our main result is that the violation of the CHSH inequality (and $n$-measurement generalisations) can be exactly characterised by a negativity witness of the hidden-variable distribution defined over the local states, and that there exists quasiprobability distributions which can saturate (up to the no-signalling limit) any such violation, whilst still having well-defined local statistics. This shows that it is possible to {recapture} the nonlocal features of {Bell experiments} through having a finite amount of negativity allowed in a hidden-variable distribution over scenarios which are, in themselves, entirely local and classical.
	\begin{SCfigure*}
	\centering
		\includegraphics[width=0.6\textwidth]{Bell_test_crop.pdf}
	\caption{A source $S$ distributes a system between two spatially separated observers, Alice ($A$) and Bob ($B$). Alice and Bob choose to measure their local part of the system with measurements $x_A,x_B$, with possible outcomes  {$y_{A},y_B\in\{-1,+1\}$}. The statistics of their measurement outcomes depending upon the physical system being distributed by the source.     \label{fig:box1}}
	\end{SCfigure*}
		\section{Setup}
	{Let us consider the following experimental setup. A source $S$ distributes a system between $2$ observers, the $k$th observer can choose some measurement $x_k\in\{0,1,\dots,X_k\}$ and record some outcome $y_k\in\{1, 2,\dots,Y_k\}$, the possible values of $k$ being $\{A,B\}$. A specific experimental {setup} is characterised by the conditional probability,
		\begin{align}
		P\left(y_A,y_B|x_A,x_B\right).
		\end{align}
		{The physical theory} governing the behaviour of the system and experiment determines the achievability of certain conditional probability distributions resulting from these experiments. We are interested in the following three physical theories:}

	$\bullet$ \textit{Classical theory} admits probability distributions of the following form,
	\begin{align}
	P\left(y_A,y_B|x_A,x_B\right)=\sum_{\lambda_A,\lambda_B}P_A\left(y_A|x_A,\lambda_A\right)P_B\left(y_B|x_B,\lambda_B\right)P_\Lambda\left(\lambda_A,\lambda_B\right),\label{eq:classical}
	\end{align}
	where $P_\Lambda\left(\lambda_A,\lambda_B\right)$ is {a} joint probability distribution defined over local hidden variables. {With each choice of hidden variables we associate a local scenario governed by ordinary local probability distributions $P_k\left(y_k|x_k,\lambda_k\right)$ for the observables $y_A,y_B$. The hidden-variable probability distribution $P_\Lambda$ determines how such local scenarios are mixed, and the probability distributions $P_k\left(y_k|x_k,\lambda_k\right)$ are called ``$\lambda_k$-local'' because they belong to the scenario associated with a particular value of $\lambda_k$, not to be confused with the (observable) marginal probability distributions that are obtained by marginalising the total probability distribution, equation \eqref{eq:classical}.}

	{$\bullet$ \textit{Quantum theory} endows us with a Hilbert space structure for our quantum states that  admits probability distributions of the following form, \begin{align}
		P\left(y_A,y_B|x_A,x_B\right)=\Tr\left[\left(M^{(A)}_{y_A|x_A}\ox M^{(B)}_{y_B|x_B}\right)\rho\right],\label{eq:prob_quantum}
		\end{align}
		where $\rho\succeq 0$ and $M^{(k)}_{y_k|x_k}$ are Positive Operator-Valued Measures (POVMs) \cite{2000quantum} for each $k$. }

	{$\bullet$ \textit{No-signalling principle}, our third physical theory, prohibits the sending of information faster than the speed of light \cite{popescu1994quantum,peres2004quantum}. Such a theory has the conditions on its conditional probability distribution that for any {$k\in \{A,B\}$},
		\begin{align}
		\sum_{y_k} P\left(y_A,y_B|x_A,x_B\right)\label{eq:cond_prob}
		\end{align} is independent of $x_k$.}
	These three physical theories ranging from the most restrictive(classical), to the least restrictive (no-signalling), with quantum theory existing somewhere between the two \cite{popescu1994quantum}.

	{Representing the full set of correlations that the quantum conditional probability distribution in equation \eqref{eq:prob_quantum} allows one to reach is a notorious problem, and the set has recently been shown to be not closed \cite{slofstra2019set}. Therefore we instead restrict ourselves to studying the achievable bounds that these conditional probabilities allow {one} to reach in nonlocal experiments {; the} original and most famous of which being the Bell inequality \cite{bell1964einstein}.
		\begin{definition}[Bell inequality] Given observers {$A$ and $B$}, each with measurement choice $x_k\in\{0_k,1_k\}$ with outcomes $y_k\in\{-1,+1\}$, experiments performed on the systems adhere to the bound,
			\begin{align}
			\left|E(0_A,0_B)-E(0_A,1_B)+E(1_A,0_B)+E(1_A,1_B)\right|\leq X,\label{equ:all_bell}
			\end{align}
			where both the correlation measure $E(x_A,x_B)=\sum_{y_A,y_B}y_Ay_B\, P(y_A,y_B|x_A,x_B)$ and bound $X\in\mathbb{R}^+$ are theory {dependent}. {The left-hand side of that inequality is often called the score of the experiment.}
			\label{def:Bell}
	\end{definition}}
	Each physical theory admits a different conditional probability distribution, and hence a different achievable bound $X$.
	{Classical theory has the CHSH bound $X=2$ \cite{clauser1969proposed}, quantum theory has the Tsirelson bound of $X=2\sqrt{2}$ \cite{cirel1980quantum} and non-signalling distributions $X=4$ \cite{popescu1994quantum}.} We are interested in the achievable bounds of a classical system's probability distribution when the hidden-variable distribution in said probability distribution can be negative.

	\section{Results}

	We now define an important object for this work, the \textit{quasiprobability distribution}.
	\begin{definition}[Quasiprobability distribution]\label{def:quasi_prob}
		We define a quasiprobability distribution as {$\tl{P}_\Lambda:\Lambda_1\times\dots\times\Lambda_N\rightarrow\mathbb{R}$} where {$\Lambda_i\subset \mathbb{R}$ and } $|\Lambda_i|<\infty\,\forall \,i$, that is properly normalised, such that,
		\begin{align}
		\sum_{\lambda_1,\dots,\lambda_N}\tl{P}_\Lambda\left(\lambda_1,\dots,\lambda_N\right)=1.
		\end{align}
	\end{definition}
	It can be seen that the collection of functions adhering to the above definition forms a convex set, which we will denote $\tl{\P}$, a super set of the convex set of positive probability distributions $\P\subset\tl{\P}$.

	We must now determine how to quantify the presence of negativity in our {quasiprobability distributions}. To this end we will use a well{-}known method for quantitatively detecting properties of a quantum state, witnesses \cite{terhal2000bell,lewenstein2000optimization,eisert2007quantitative,chabaud2021witnessing}. Let us therefore proceed by defining a \emph{negativity witness:}
	\begin{definition}[Negativity witness]\label{def:neg_witness} Given some properly normalised probability distribution $P$, a well-defined \textit{negativity witness} is one which,
		\begin{align}
		\mathcal{N}(P)=0\,\,\forall\,P\in{\P}.
		\end{align}
		We may additionally require such a witness to `faithfully' detect negativity,
		\begin{align}
		\mathcal{N}(P)>0\,\,\forall\,P\in \tilde{\P}{\backslash}\P.
		\end{align}
	\end{definition}

	In the following, we consider classical local hidden-variable models{ as defined in } equation \eqref{eq:classical}, but we replace the hidden-variable probability distribution $P_\Lambda$ with a quasiprobability distribution $\tilde{P}_\Lambda$,
		\begin{align}
		P\left(y_A,y_B|x_A,x_B\right)=\sum_{\lambda_A,\lambda_B}P_1\left(y_A|x_A,\lambda_A\right)P_2\left(y_B|x_B,\lambda_B\right)\tilde{P}_\Lambda\left(\lambda_A,\lambda_B\right).\label{eq:quasi_hidden}
		\end{align}

		This corresponds to a scenario where different local statistics of observations, governed by {$\lambda$-local (ordinary)} probability distributions $P_k\left(y_k|x_k,\lambda_k\right)$, are mixed according to a quasiprobability distribution $\tilde{P}_\Lambda$.
		{However, when $\tilde{P}_\Lambda$ takes negative values, we should no longer think of the model as an ignorance mixture of valid local scenarios but rather as a nonlocal model \cite{al2013simulating}. Furthermore, when compared with ordinary hidden-variable models not all combinations of hidden-variable and $\lambda$-local probability distributions are valid; only those combinations  which lead to well-defined $	P\left(y_A,y_B|x_A,x_B\right)$ {are valid}, i.e., comprised of values between 0 and 1 (the normalisation condition is always fulfilled).}

 In addition to the correlation function  between {two measurements $x_A$ and $x_B$,
${E}(x_A,x_B)$, it will also be useful to define $\lambda_k$-local expectation values {corresponding to an imagined scenario where} observer $k$ {is able to perform} measurement $x_k$ in {the} local scenario {corresponding to} $\lambda_k$,
		\begin{align}
		\inner{k}_{\lambda_k}^{x_k}:=\sum_{y_k}y_k P_k(y_k|x_k,\lambda_k). \label{eq:local_exp}
		\end{align}
		This $\lambda_k$-local expectation value will be useful to formulate our results, but does not correspond to the actual observations which are themselves governed by equation \eqref{eq:quasi_hidden}. }

	We are now in a position to state the main result of this work{, the quasiprobabilistic Bell inequality}.
	\begin{theorem}[{Quasiprobabilistic Bell inequality}]\label{thm:CHSH_break}
		Given observers {$A$ and $B$}, each with measurement choice $x_k\in\{0_k,1_k\}$ with outcomes $y_k\in\{-1,+1\}$ whose systems are distributed according to some  {quasi}probability  distribution $\tilde{P}_\Lambda$, then the quasiprobabilistic Bell inequality holds:
		\begin{align}
		\left|{E}(0_A,0_B)-{E}(0_A,1_B)+{E}(1_A,0_B)+{E}(1_A,1_B)\right|\leq\,2+\mathcal{N}(\tilde{P}_\Lambda),\label{equ:quasi_CHSH2}
		\end{align}
		\color{black}
		{where}
		{\begin{align}
			\mathcal{N}(\tilde{P}_\Lambda)\coloneq\begin{cases}
			\mathcal{N}_+(\tilde{P}_\Lambda) \qquad\text{ if }{E}(1_A,0_B)+{E}(1_A,1_B)<0,\\
			\mathcal{N}_-(\tilde{P}_\Lambda) \qquad\text{ else,}
			\end{cases}
			\end{align}
		{is a negativity witness, and }
			 $$\mathcal{N}_\pm(\tilde{P}_\Lambda):=\sum_{\lambda_A,\lambda_B}\left[2\pm\left(\inner{A}_{\lambda_A}^{1_A}\inner{B}_{\lambda_B}^{1_B}+\inner{A}_{\lambda_A}^{1_A}\inner{B}_{\lambda_B}^{0_B}\right)\right]\left(\Big{|}\tilde{P}_\Lambda(\lambda_A,\lambda_B)\Big{|}-\tilde{P}_\Lambda(\lambda_A,\lambda_B) \right).$$  }
	\end{theorem}
	The proof of this theorem begins analogously with Bell's proof of the CHSH bound \cite{bell2004speakable}, but diverges when the assumption $P\in\P$ is {made in Bell's proof}. The above result shows that if an arbitrary amount of negativity is allowed in  {the hidden-variable probability distribution} then the upper bound of equation \eqref{equ:quasi_CHSH2}  can be arbitrarily large. However, it should be noted that a natural limit of $4$ in the relevant Bell tests (i.e., for the upper bound in the quasiprobabilistic Bell inequality) is imposed by the requirement that $	P\left(y_A,y_B|x_A,x_B\right)$ is a well-defined, valid probability distribution \footnote{This can be easily seen as for {well-defined experimental setups where $y_k=\{-1,+1\}$, i.e.~for valid} $	P\left(y_A,y_B|x_A,x_B\right)$, each term on the left hand side of equation \eqref{equ:quasi_CHSH2} can be at most $1$ and least $-1$. }.

	The previous result of Al-Safi and Short \cite{al2013simulating} {showed} that it was possible to violate said inequality up to this no-signalling bound of $X=4$. Therefore, in order to emulate the physical results seen in Bell test{s} (Tsirelson bound) one needs a negative probability distribution whose witness equals $\mathcal{N}(\tilde{P}_\Lambda)=2(\sqrt{2}-1)$. In section \ref{sec:examples} we show that for any $\mathcal{N}({\tilde{P}_\Lambda}){\leq 2}$, there exist {quasiprobabilistic hidden-variable models with valid local measurement statistics that}  saturate inequality \eqref{equ:quasi_CHSH2}. We would hope that if a physical mechanism was discovered that allowed a joint hidden-variable probability distribution to have the appearance of negativity, one would expect that said physical mechanism was limited in such a way that it resulted in the Tsirelson bound and more generally was able to reconstruct the limits on quantum correlations.

	It is also important to note that although said witness $\mathcal{N}(\tilde{P}_\Lambda)$ is a valid witness according to definition \ref{def:neg_witness} it is not necessarily a `faithful' one. However this can be rectified, at the cost of loosening the bound, by redefining said witness. For example the function {$\mathcal{N}'({\tilde{P}_\Lambda}):=\sum_{\lambda_A,\lambda_B}4\left(\Big{|}\tilde{P}_\Lambda(\lambda_A,\lambda_B)\Big{|}-\tilde{P}_\Lambda(\lambda_A,\lambda_B) \right)$}, is defined to be both a valid and `faithful' witness.

	There are numerous generalisations of the famous CHSH inequalities, {such as} multiple parties \cite{perez2008unbounded}, arbitrary outcomes \cite{cope2019bell}, etc. These would no doubt be interesting to study but we leave it to future work to explore these other generalisations and instead focus on the scenario in which Alice and Bob have access to an arbitrary number of measurement settings \cite{wehner2006tsirelson}.
	\begin{theorem}\label{thm:CHSH_break_N}
		Given observers {$A$ and $B$}, each with $n\geq2$ measurements
		$x_k\in\{0_k,1_k,\dots,{n-1}_k\}$ with outcomes $y_k\in\{-1,+1\}$ whose systems are distributed according to some  {quasi}probability distribution $\tilde{P}_\Lambda$,
		{\begin{align}
			\left|\sum^{n-1}_{i=0}{E}(i_A,i_B)+\sum^{n-1}_{i=1}{E}(i_A,i-1_B)-{E}(0_A,n-1_B)\right|\leq\,2n-2+{\mathcal{N}_n({\tilde{P}_\Lambda})}, \label{eq:CHSH_break_N}
			\end{align}}
		where { $\mathcal{N}_n({\tilde{P}_\Lambda})=\sum^{n-1}_{i=1}\mathcal{N}^{(i)}(\tilde{P}_\Lambda)$ {is a negativity witness} with}
		\begin{align}
		\mathcal{N}^{(x)}(\tilde{P}_\Lambda)\coloneq\begin{cases}
		\mathcal{N}^{(x)}_+(\tilde{P}_\Lambda) \qquad\text{ if }{E}(0_A,x_B)+{E}(0_A,x-1_B)<0,\\
		\mathcal{N}^{(x)}_-(\tilde{P}_\Lambda) \qquad\text{ else,}
		\end{cases}\label{eq:neg_cases_n}
		\end{align}
		where $\mathcal{N}^{(x)}_{\pm}({\tilde{P}_\Lambda})\coloneqq\sum_{\lambda_A,\lambda_B}\left[2\pm\left(\inner{A}_{\lambda_A}^{x_A}\inner{B}_{\lambda_B}^{x_B}+\inner{A}_{\lambda_A}^{x_A}\inner{B}_{\lambda_B}^{x-1_B}\right)\right]\left(\Big{|}{\tilde{P}_\Lambda}(\lambda_A,\lambda_B)\Big{|}-{\tilde{P}_\Lambda}(\lambda_A,\lambda_B) \right)$.
	\end{theorem}
	The proof of the above theorem can be found in appendix \ref{app:proof_thm2}, it utilises proof by induction by chaining together the inequalities from theorem \ref{thm:CHSH_break}. In section \ref{sec:examples} we show that the bound in theorem \ref{thm:CHSH_break_N} can be saturated. Namely, {for any $\mathcal{N}_n(\tilde{P}_\Lambda){\leq 2}$, there exist well-defined $	P\left(y_A,y_B|x_A,x_B\right)$, characterised by a  quasiprobability hidden-variable distribution $\tilde{P}_\Lambda(\lambda_A,\dotsc,\lambda_n)$, that saturate inequality \eqref{eq:CHSH_break_N}.} In addition, analogously to the two measurement result, at the cost of loosening the bound we can ensure that the above witness is also `faithful' by  choosing for all $x$, $\mathcal{N}^{(x)}(\tilde{P}_\Lambda)=\mathcal{N}'({\tilde{P}_\Lambda})$.

	\section{Example}\label{sec:examples}
	In order to understand how to saturate the Bell inequality from Theorem \ref{thm:CHSH_break}, we rewrite the left-hand side of equation \eqref{equ:quasi_CHSH2} as,
	\begin{align}
\Bigg{|}\sum_{\lambda}\M(\lambda)\,\tilde{P}_\Lambda(\lambda)\Bigg{|}. \label{eq:CHSH_example2}
	\end{align}
We replaced the hidden variables $\lambda_A$ and $\lambda_B$ with a single hidden variable $\lambda$ because our example only uses a single hidden variable $\lambda$. Further,   $\M (\lambda):=\inner{A}_{\lambda}^{0_A}\inner{B}_{\lambda}^{0_B}-\inner{A}_{\lambda}^{0_A}\inner{B}_{\lambda}^{1_B}+\inner{A}_{\lambda}^{1_A}\inner{B}_{\lambda}^{0_B}+\inner{A}_{\lambda}^{1_A}\inner{B}_{\lambda}^{1_B}$ are the scores of each of the $\lambda$-local distributions {, that is }  $-2\leq \M(\lambda)\leq 2$ {holds}.

	For a given value of the negativity witness, we exceed the local bound maximally by the simple strategy of weighting classical distributions with $\M (\lambda)=+2$ with positive quasiprobability, while simultaneously taking a classical distribution with $\M (\lambda)=-2$ with {negative weight}.
    To ensure that the {total}  probability distribution $P\left(y_A,y_B|x_A,x_B\right)$ is well-defined we {make} a choice of three deterministic classical distributions with positive weight and a fourth with negative weight. Our four deterministic classical distributions can be denoted $[(-,-)_A,(-,+)_B], [(+,-)_A,(-,-)_B], [(+,+)_A,(+,+)_B]$ and $[(+,-)_A,(-,+)_B]$. { Here, our notation means} that the distributions can be produced by assigning the first pair of symbols to Alice and the second to Bob. Each party chooses to read either the first or second of the symbols given to them (this choice reflects their measurement setting $x_k$) while the outcome of their measurement is determined by the symbol itself; that is{,} $y_k=+1$ ($y_k=-1$) for a plus (minus) sign. This {experimental description of distributing classical information makes clear that these distributions are local}{,} with our hidden variable $\lambda$ indicating which of these sets the source actually produces.

    The source produces each of the distributions according to the following quasiprobability distribution,
         \begin{align}\label{eq:prob_sat_2}
         \tilde{P}_\Lambda(\lambda)=\begin{cases}
 	\frac{4+\mathcal{N}}{12}\qquad\text{ for }\lambda=1,2,3,\\
 	-\frac{\mathcal{N}}{4} \qquad\,\,\,\text{  for }\lambda=4,
 		\end{cases}
     \end{align}
    where{ $\M (\lambda)=2$ if $\lambda=1,2,3$  and  $\M (\lambda)=-2$ if $\lambda=4$}.  We can use tables to represent $\lambda$-local probability distributions, and the total probability distribution is then given as the weighted sum of such tables:

\begin{align}\label{CHSH_table}
&\scriptsize	\frac{4 + \mathcal{N}}{12} \begin{pmatrix}                    \begin{bmatrix}[c|cccc]
		     & \multicolumn{4}{c}{y_A y_B} \\
			x_A x_B & -- & -+ & +- & ++ \\ \hline
			00        & 1 & 0 & 0 & 0 \\
			01        & 0 & 1 & 0 & 0 \\
			10        & 1 & 0 & 0 & 0 \\
			11        & 0 & 1 & 0 & 0 \\
		\end{bmatrix}&+\begin{bmatrix}[c|cccc]
		     & \multicolumn{4}{c}{y_A y_B} \\
			x_A x_B & -- & -+ & +- & ++ \\ \hline
			00        & 0 & 0 & 1 & 0 \\
			01        & 0 & 0 & 1 & 0 \\
			10        & 1 & 0 & 0 & 0 \\
			11        & 1 & 0 & 0 & 0 \\
		\end{bmatrix}+\begin{bmatrix}[c|cccc]
		     & \multicolumn{4}{c}{y_A y_B} \\
			x_A x_B & -- & -+ & +- & ++ \\ \hline
			00        & 0 & 0 & 0 & 1 \\
			01        & 0 & 0 & 0 & 1 \\
			10        & 0 & 0 & 0 & 1 \\
			11        & 0 & 0 & 0 & 1 \\
		\end{bmatrix}
	 \end{pmatrix}
	\scriptsize	- \frac{\mathcal{N}}{4}			\begin{bmatrix}[c|cccc]
		     & \multicolumn{4}{c}{y_A y_B} \\
			x_A x_B & -- & -+ & +- & ++ \\ \hline
			00        & 0 & 0 & 1 & 0 \\
			01        & 0 & 0 & 0 & 1 \\
			10        & 1 & 0 & 0 & 0 \\
			11        & 0 & 1 & 0 & 0 \\
		\end{bmatrix}
		\nonumber\\ 	\nonumber\\
	 &\scriptsize\qquad\quad\,\,	=\frac{1}{12}\begin{bmatrix}[c|cccc]
		     & \multicolumn{4}{c}{y_A y_B} \\
			x_A x_B & -- & -+ & +- & ++ \\ \hline
	00 & 4+\mathcal{N} & 0 & 4-2\mathcal{N} & 4+\mathcal{N} \\
	01 &	     0 & 4+\mathcal{N} & 4+\mathcal{N} &  4-2\mathcal{N} \\
	10 &	     8 - \mathcal{N} & 0 & 0 & 4+\mathcal{N} \\
	11 &		 4+\mathcal{N} & 4-2\mathcal{N} & 0 & 4+\mathcal{N} \\
		\end{bmatrix}.
	\end{align}
The requirement that the resulting total probability distribution must be valid implies
 $\mathcal{N}\leq2$ {which corresponds to the no-signalling limit.} {Furthermore, it is easy to check that said distribution indeed gives a value of $\mathcal{N}$ for the negativity witness.}

The quasiprobabilistic Bell inequality score for this experiment is $ 2+\mathcal{N}$, which upon substituting equation \eqref{eq:prob_sat_2} into the negativity witness, can be seen to saturate the bound. In appendix \ref{app:sat_n_ex} we discuss how one can generalise the above to the $n$-measurement scenario.

	\section{Conclusion}
	We have shown that there exists a relationship between the amount of negativity allowed in a joint hidden-variable distribution, and the degree to which said distribution can demonstrate nonlocality in a Bell experiment. In particular, theorem  \ref{thm:CHSH_break_N} {introduces a quasiprobabilistic Bell inequality,} which gives us a sharp bound in the scenario of two parties with {$n$} inputs (corresponding to a choice between $n$ measurements) and can be used straightforwardly to reconstruct quantum statistics using nothing more than local, separable classical probability distributions and a quasiprobability distribution over them -- granted an appropriately well spent budget of negativity.

	Our work sits within the long-established tradition of trying to understand quantum theory through interpretative lenses which remove some particular aspect from a classical worldview. Such approaches are wide and varied, including superdeterminism \cite{hossenfelder2020rethinking, adlam2018quantum}; retro-causality \cite{price2012does}; invoking an irreducible role for subjectivity in physics \cite{mermin2014physics, fuchs2014introduction, mueller2020law}; taking physical reality to consist of interacting, separate realms \cite{goldstein2001bohmian, esfeld2014ontology}; allowing the relativity of pre and post-selection \cite{bacciagaluppi2020reverse}; taking Hilbert space to be literal \cite{carroll2021reality}, and so on. Here we add to this list, in that we present an additional way to re-capture the nonlocal features of quantum theory: through having a finite amount of negativity allowed in a hidden-variable distribution over scenarios which are, in themselves, entirely local and classical.  We are not claiming that such quasi distributions are `real' -- only, more modestly, that such a perspective could not be ruled out at this stage.

	Pursuing this line of reasoning, we would hope that our results may help to determine the fundamental restrictions on a system's quasiprobability hidden-variable distribution such that it captures the full character of  physical correlations. Put another way; we know that zero negativity can capture the set of classical correlations, whilst un-bounded negativity can capture the non-signalling set. Given that the set of quantum correlations lies between these two -- what are the restrictions on the quasiprobability hidden-variable distribution which would suffice to identify the full set of quantum correlations? We hope to explore this question in further work.

	\acknowledgements{We are grateful to Benjamin Yadin, Richard Moles and Thomas Veness for helpful discussions and the Physical Institute for Theoretical Hierarchy (PITH) for encouraging an investigation into this topic. BM acknowledges financial support from the Engineering and Physical Sciences Research Council (EPSRC) under the Doctoral Prize Grant (Grant No.~EP/T517902/1). LF acknowledges financial support from the Austrian Science Fund (FWF) through SFB BeyondC (Grant No.~F7102). BL is supported by Leverhulme Trust Research Project Grant (RPG-2018-213). D.G. acknowledges support from FQXI (grant no. RFP-IPW-1907)}