\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.
  
