# A new coupled three-form dark energy model and implications for the $H_{0}$
tension

Yan-Hong Yao<EMAIL_ADDRESS>Xin-He Meng<EMAIL_ADDRESS>Department
of Physics, Nankai University, Tianjin 300071, China

###### Abstract

We propose a new coupled three-form dark energy model to relieve the Hubble
tension in this paper. Firstly, by performing a dynamical analysis with the
coupled three-form dark energy model, we obtain four fixed points, including a
saddle point representing a radiation dominated Universe, a saddle point
representing a matter dominated Universe, and two attractors representing two
saturated de Sitter Universes. Secondly, by confronting the coupled three-form
dark energy model and the $\Lambda$ cold dark matter model (the $\Lambda$CDM
model) with cosmic microwave background (CMB), baryonic acoustic oscillations
(BAO), Type Ia supernovae (SN Ia) observations, we obtain
$H_{0}=67.8_{-0.6}^{+0.7}$($1\sigma$ level) km/s/Mpc for the coupled three-
form dark energy model and $H_{0}=67.6_{-0.5}^{+0.5}$($1\sigma$ level)
km/s/Mpc for the $\Lambda$CDM model, the former is in strong tension with the
latest local measured $H_{0}$ value at $4.3\sigma$ confidence level, while the
latter is in strong tension with the latest local measured $H_{0}$ value at
$5.1\sigma$ level.

##  I Introduction

The dominant components of the Universe, dark energy and dark matter, are
incorporated in the standard model of cosmology, known as the $\Lambda$CDM
model. In such model, dark energy takes the form of a cosmological constant
$\Lambda$. Meanwhile, the dark matter is nonrelativistic, and it interacts
with ordinary matter only through gravity.

Although mathematically simple, the $\Lambda$CDM model provides an excellent
fit to a wide range of cosmological data. However, an exception is emerging in
the Hubble constant $H_{0}$. In 2018, the Planck satellite measured a $H_{0}$
value of $67.4\pm 0.5$($1\sigma$ level) km/s/Mpc from a $\Lambda$CDM fit to
the CMB aghanim2018planck . In 2019, The SH0ES collaboration yielded a latest
$H_{0}$ value of $74\pm 1.4$($1\sigma$ level) km/s/Mpc from direct
measurements by using so-called standard candles: type Ia supernovae and
Cepheid variable stars riess2019large . Recently, the H0LiCOW collaboration
obtained a $H_{0}$ value of $73.3_{-1.8}^{+1.7}$($1\sigma$ level) km/s/Mpc by
using gravitationally lensed quasar wong2019h0licow . Combining the SH0ES and
H0LiCOW measurements gives a model independent $H_{0}$ value of $73.8\pm
1.1$($1\sigma$ level) km/s/Mpc, which is in 5.3$\sigma$ tension with the
$\Lambda$CDM prediction.

Since preliminary attempt to resolve the $H_{0}$ tension by searching
considerable systematic errors in the Planck observation and the local
measurements have failed spergel2015planck ; addison2016quantifying ;
aghanim2017planck ; cardona2017determining ; follin2018insensitivity ,
increasing attention is focusing on the possibility that the $\Lambda$CDM
model is not the final picture. For example, ref.Li2013Planck ; huang2016how
pointed out that a phantom dark energy prefers a high value of $H_{0}$. Refs.
battye2014evidence ; Zhang2014Neutrinos ; zhang2015sterile ; feng2018searching
; zhao2018measuring ; choudhury2019constraining  showed that considering extra
relativistic degrees of freedom $N_{eff}$ in the $\Lambda$CDM model favors a
high value of $H_{0}$ when $N_{eff}>3.046$. Although phantom dark energy and
extra relativistic species can help with the $H_{0}$ tension, it is worth to
mention that these solutions are disfavored from both BAO and SN Ia data and
from a model comparison point of view. vagnozzi2020new  Refs. agrawal2019rock
; poulin2019early  introduced an exotic early dark energy (EDE) that acts as a
cosmological constant before a critical redshift $z_{c}$ around 3000 but whose
density then dilutes faster than radiation to resolve the $H_{0}$ tension.
Ref.li2019simple  introduced a emergent dark energy with its equation of state
increases from $-\tfrac{2}{3\mathrm{ln}\,10}-1$ in the past to $-1$ in the
future to handle the Hubble problem. Ref.2018Vacuum  showed that the Parker
vacuum metamorphosis (VM) model, physically motivated by quantum gravitational
effects, can remove the Hubble tension. In addition, an interacting dark
energy can affect the constraint results of $H_{0}$, which provides another
way to address the Hubble constant problem.di2017can ; yang2018interacting

In Ref.yao2018a , we proposed a power-law coupled three-form dark energy model
and successfully alleviated the coincidence problem with this model. Agreeing
with the opinion in yao2018a  that dark energy might be represented by a
three-form field, and considering the fact that interacting dark energy
affects the constraint results of $H_{0}$, in this paper, we put forward a new
coupled three-form dark energy model to relieve the Hubble tension.

The contents of this paper are as follows. In section II, we consider a
coupled dark energy model in which dark energy is represented by a three-form
field and other components are represented by ideal fluids. In section III, we
derive the autonomous system of evolution equations, and analyze the stability
of its fixed points. In section IV, we confront the model with the data from
CMB, BAO, SN Ia observations. In the last section, we make a brief conclusion
with this paper. For convenience, we set 100 km/s/Mpc=1,i.e., $H_{0}=h$ in the
following part of the paper.

##  II A coupled three-form dark energy model

In this section, a new coupled three-form dark energy model is presented, in
which dark energy is represented by a three-form field and other components
are represented by ideal fluids. We restrict the coupling here to be the
conformal form koivisto2013coupled ; yao2018a , a case that has been
thoroughly studied in the context of scalar fields. The total Lagrangian is
written as

$\mathcal{L}=\frac{R}{2\kappa^{2}}+\mathcal{L}_{(d)}+\mathcal{L}_{(b)}+\mathcal{L}_{(\gamma)}+\mathcal{L}_{(\nu)},$
(1)

where $R$ denotes the Ricci scalar and $\kappa=\sqrt{8\pi G}$ is the inverse
of the reduced Planck mass. $d$, $b$, $\gamma$, and $\nu$ denote dark sectors,
baryon, photon, and neutrino, respectively. For ideal fluids, the derivation
of the equations of motion from a variational principle is complicated since
one needs to consider the constraint equations satisfied by the fluid
variables. To solve this problem, several variational formulations have been
proposed. In this paper, we consider the variational formulations discussed in
ray1972lagrangian , then each Lagrangian can be expressed as

$\displaystyle\mathcal{L}_{(d)}$ $\displaystyle=$
$\displaystyle-\frac{1}{48}F^{2}-I(A^{2})\tilde{\rho}_{(c)}+\lambda_{1}(g_{\mu\nu}u_{(c)}^{\mu}u_{(c)}^{\nu}+1)+\lambda_{2}\nabla_{\alpha}(\tilde{\rho}_{(c)}u_{(c)}^{\alpha}),$
(2) $\displaystyle\mathcal{L}_{(b)}$ $\displaystyle=$
$\displaystyle-\rho_{(b)}+\lambda_{3}(g_{\mu\nu}u_{(b)}^{\mu}u_{(b)}^{\nu}+1)+\lambda_{4}\nabla_{\alpha}(\rho_{(b)}u_{(b)}^{\alpha}),$
(3) $\displaystyle\mathcal{L}_{(\gamma)}$ $\displaystyle=$
$\displaystyle-\tilde{\rho}_{(\gamma)}(1+\epsilon_{(\gamma)}(\tilde{\rho}_{(\gamma)}))+\lambda_{5}(g_{\mu\nu}u_{(\gamma)}^{\mu}u_{(\gamma)}^{\nu}+1)+\lambda_{6}\nabla_{\alpha}(\tilde{\rho}_{(\gamma)}u_{(\gamma)}^{\alpha}),$
(4) $\displaystyle\mathcal{L}_{(\nu)}$ $\displaystyle=$
$\displaystyle-\tilde{\rho}_{(\nu)}(1+\epsilon_{(\nu)}(\tilde{\rho}_{(\nu)}))+\lambda_{7}(g_{\mu\nu}u_{(\nu)}^{\mu}u_{(\nu)}^{\nu}+1)+\lambda_{8}\nabla_{\alpha}(\tilde{\rho}_{(\nu)}u_{(\nu)}^{\alpha}),$
(5)

$A$ and $F=dA$ represent the three-form field and the field strength tensor,
$I(A^{2})$ is referred to as coupling function. $\tilde{\rho}$ denotes the
rest density, and $\epsilon$ denotes the rest, specific internal energy, which
is connected with pressure through the following relation,

$\displaystyle\frac{d\epsilon}{d\tilde{\rho}}$ $\displaystyle=$
$\displaystyle\frac{p}{\tilde{\rho}^{2}}$ (6)

Finally, $\lambda_{1}-\lambda_{8}$ are the multipliers.

One now can obtain the field equations from the total action

$S=\int\mathcal{L}\sqrt{-g}d^{4}x.$ (7)

The variation of the action with respect to the
$g^{\mu\nu},\tilde{\rho}_{(c)},u_{(c)}^{\alpha},\rho_{(b)},u_{(b)}^{\alpha},\tilde{\rho}_{(\gamma)},u_{(\gamma)}^{\alpha},\tilde{\rho}_{(\nu)},u_{(\nu)}^{\alpha}$
leads to

$\displaystyle R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R$ $\displaystyle=$
$\displaystyle\kappa^{2}T_{\mu\nu},$ (8)
$\displaystyle(\partial_{\alpha}\lambda_{2})u_{(c)}^{\alpha}$ $\displaystyle=$
$\displaystyle-I,$ (9) $\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}I\tilde{\rho}_{(c)},$ (10)
$\displaystyle(\partial_{\alpha}\lambda_{4})u_{(b)}^{\alpha}$ $\displaystyle=$
$\displaystyle-1,$ (11) $\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\rho_{(b)},$ (12)
$\displaystyle(\partial_{\alpha}\lambda_{6})u_{(\gamma)}^{\alpha}$
$\displaystyle=$
$\displaystyle-(1+\epsilon_{(\gamma)}(\tilde{\rho}_{(\gamma)})),$ (13)
$\displaystyle\lambda_{5}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\epsilon_{(\gamma)}(\tilde{\rho}_{(\gamma)}))\tilde{\rho}_{(\gamma)},$
(14) $\displaystyle(\partial_{\alpha}\lambda_{8})u_{(\nu)}^{\alpha}$
$\displaystyle=$ $\displaystyle-(1+\epsilon_{(\nu)}(\tilde{\rho}_{(\nu)})),$
(15) $\displaystyle\lambda_{7}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\epsilon_{(\nu)}(\tilde{\rho}_{(\nu)}))\tilde{\rho}_{(\nu)},$
(16)

with the help of equation (6), and $(9)-(16)$, the total energy-momentum
tensor for three-form field and ideal fluids is written as

$\begin{split}T_{\mu\nu}=&\frac{1}{6}F_{\mu\alpha\beta\gamma}F_{\nu}^{\alpha\beta\gamma}+6\frac{dlnI}{dA^{2}}\rho_{(c)}A_{\mu}^{\alpha\beta}A_{\nu\alpha\beta}-g_{\mu\nu}\frac{1}{48}F^{2}\\\
&+\rho_{(c)}u_{(c)\mu}u_{(c)\nu}+\rho_{(b)}u_{(b)\mu}u_{(b)\nu}+(\rho_{(\gamma)}+p_{(\gamma)})u_{(\gamma)\mu}u_{(\gamma)\nu}+p_{(\gamma)}g_{\mu\nu}+(\rho_{(\nu)}+p_{(\nu)})u_{(\nu)\mu}u_{(\nu)\nu}+p_{(\nu)}g_{\mu\nu},\end{split}$
(17)

where $\rho_{(c)}=I(A^{2})\tilde{\rho}_{(c)}$,
$\rho_{(\gamma)}=\tilde{\rho}_{(\gamma)}(1+\epsilon_{(\gamma)}(\tilde{\rho}_{(\gamma)}))$,
and
$\rho_{(\nu)}=\tilde{\rho}_{(\nu)}(1+\epsilon_{(\nu)}(\tilde{\rho}_{(\nu)}))$.
Since we assume that the Universe is homogeneous and isotropic, we have
$u_{\mu}=u_{(c)\mu}=u_{(b)\mu}=u_{(\gamma)\mu}=u_{(\nu)\mu}$. In addition, we
denote $\rho_{(r)}=\rho_{(\gamma)}+\rho_{(\nu)}$ and
$p_{(r)}=p_{(\gamma)}+p_{(\nu)}$. As a result, (17) becomes to

$\begin{split}T_{\mu\nu}=&\frac{1}{6}F_{\mu\alpha\beta\gamma}F_{\nu}^{\alpha\beta\gamma}+6\frac{dlnI}{dA^{2}}\rho_{(c)}A_{\mu}^{\alpha\beta}A_{\nu\alpha\beta}-g_{\mu\nu}\frac{1}{48}F^{2}\\\
&+\rho_{(c)}u_{\mu}u_{\nu}+\rho_{(b)}u_{\mu}u_{\nu}+(\rho_{(r)}+p_{(r)})u_{\mu}u_{\nu}+p_{(r)}g_{\mu\nu}.\end{split}$
(18)

By varying the total action with respect to the three-form field, we have the
following equation of motion,

$\nabla_{\alpha}F^{\alpha\mu\nu\rho}=12\frac{dlnI}{dA^{2}}\rho_{(c)}A^{\mu\nu\rho}.$
(19)

Using the equation of motion for the three-form field and the condition that
the divergence of the total stress energy tensor vanishes, we have the
equation of motion for the dark matter:

$\nabla_{\mu}(\rho_{(c)}u^{\mu}u_{\nu})=-2\frac{dlnI}{dA^{2}}\rho_{(c)}A^{\alpha\beta\gamma}\nabla_{\nu}A_{\alpha\beta\gamma}.$
(20)

The homogeneous, isotropic, and spatially flat space-time is described by the
Friedmann- Robertson-Walker (FRW) metric,

$ds^{2}=-dt^{2}+a(t)^{2}d\vec{x}^{2},$ (21)

and the three-form field is assumed as a time-like component of a dual vector
field in order to be compatible with FRW symmetries Koivisto2009Inflation2 .

$A_{ijk}=X(t)a(t)^{3}\varepsilon_{ijk}.$ (22)

To specify a coupled three-form dark energy model, we set the coupling
function $I$ to be

$I=(1+\frac{\kappa^{2}}{6}A^{2})^{\frac{\lambda}{2}}=(1+\kappa^{2}X^{2})^{\frac{\lambda}{2}},$
(23)

where $\lambda$ is the coupling constant. Function $I$ is selected without any
fundamental reason. Now we have the Friedmann equations:

$\displaystyle H^{2}$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}}{3}\rho,$ (24) $\displaystyle\dot{H}$
$\displaystyle=$ $\displaystyle-\frac{\kappa^{2}}{2}(\rho+p),$ (25)

where

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle-
T_{0}^{0}=-g^{00}T_{00}=\frac{1}{2}(\dot{X}+3HX)^{2}+\rho_{(c)}+\rho_{(b)}+\rho_{(r)},$
(26) $\displaystyle p$ $\displaystyle=$
$\displaystyle\frac{1}{3}T_{i}^{i}=\frac{1}{3}g^{ii}T_{ii}=-\frac{1}{2}(\dot{X}+3HX)^{2}+\lambda\frac{\kappa^{2}X^{2}}{1+\kappa^{2}X^{2}}\rho_{(c)}+\frac{1}{3}\rho_{(r)}.$
(27)

The independent equation of motion of the three-form field is

$\ddot{X}+3\dot{H}X+3H\dot{X}+\lambda\frac{\kappa^{2}X}{1+\kappa^{2}X^{2}}\rho_{(c)}=0,$
(28)

and the energy conservation equations of two dark sectors is written as

$\displaystyle\dot{\rho}_{(c)}+3H\rho_{(c)}$ $\displaystyle=$
$\displaystyle\delta H\rho_{(c)},$ (29)
$\displaystyle\dot{\rho}_{(X)}+3H(\rho_{(X)}+p_{(X)})$ $\displaystyle=$
$\displaystyle-\delta H\rho_{(c)},$ (30)

where

$\displaystyle\delta$ $\displaystyle=$
$\displaystyle\lambda\frac{\kappa^{2}XX^{\prime}}{1+\kappa^{2}X^{2}},$ (31)
$\displaystyle\rho_{(X)}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\dot{X}+3HX)^{2},$ (32) $\displaystyle p_{(X)}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}(\dot{X}+3HX)^{2}+\lambda\frac{\kappa^{2}X^{2}}{1+\kappa^{2}X^{2}}\rho_{(c)}.$
(33)

##  III Dynamical system of the coupled three-form dark energy model

In order to study the dynamical behaviors of the coupled three-form dark
energy model, it is convenient to introduce the following dimensionless
variable Koivisto2009Inflation2

$x_{1}=\kappa X,\hskip
28.45274ptx_{2}=\frac{\kappa}{\sqrt{6}}(X^{\prime}+3X),\hskip
28.45274ptx_{3}=\frac{\kappa\sqrt{\tilde{\rho}_{(c)}}}{\sqrt{3}H},\hskip
28.45274ptx_{4}=\frac{\kappa\sqrt{\rho_{(b)}}}{\sqrt{3}H}.$ (34)

The autonomous system of evolution equations then can be written as follows by
applying the Friedmann equations and equation of motion,

$\displaystyle x_{1}^{\prime}$ $\displaystyle=$
$\displaystyle\sqrt{6}x_{2}-3x_{1},$ (35) $\displaystyle x_{2}^{\prime}$
$\displaystyle=$ $\displaystyle(\frac{3}{2}(1+\frac{\lambda
x_{1}^{2}}{1+x_{1}^{2}})(1+x_{1}^{2})^{\frac{\lambda}{2}}x_{3}^{2}+\frac{3}{2}x_{4}^{2}+2(1-x_{2}^{2}-(1+x_{1}^{2})^{\frac{\lambda}{2}}x_{3}^{2}-x_{4}^{2}))x_{2}-\frac{\sqrt{6}}{2}\frac{\lambda
x_{1}}{1+x_{1}^{2}}(1+x_{1}^{2})^{\frac{\lambda}{2}}x_{3}^{2},\hskip
28.45274pt$ (36) $\displaystyle x_{3}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}x_{3}+(\frac{3}{2}(1+\frac{\lambda
x_{1}^{2}}{1+x_{1}^{2}})(1+x_{1}^{2})^{\frac{\lambda}{2}}x_{3}^{2}+\frac{3}{2}x_{4}^{2}+2(1-x_{2}^{2}-(1+x_{1}^{2})^{\frac{\lambda}{2}}x_{3}^{2}-x_{4}^{2}))x_{3},$
(37) $\displaystyle x_{4}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}x_{4}+(\frac{3}{2}(1+\frac{\lambda
x_{1}^{2}}{1+x_{1}^{2}})(1+x_{1}^{2})^{\frac{\lambda}{2}}x_{3}^{2}+\frac{3}{2}x_{4}^{2}+2(1-x_{2}^{2}-(1+x_{1}^{2})^{\frac{\lambda}{2}}x_{3}^{2}-x_{4}^{2}))x_{4}.$
(38)

Here and in the following, the prime stands for the derivative with respect to
e-folding time.

| $(x_{1},x_{2},x_{3},x_{4})$ | Existence | $\Omega_{(X)}$ | $\Omega_{(c)}$ | $\Omega_{(b)}$ | $\Omega_{(r)}$ | $w_{(X)}$ | $w_{eff}$
---|---|---|---|---|---|---|---|---
(a) | $(0,0,0,0)$ | All | $0$ | $0$ | 0 | $1$ | $-$ | $\frac{1}{3}$
(b) | $(0,0,x_{3},\sqrt{1-x_{3}^{2}})$ | $0\leq x_{3}\leq 1$ | $0$ | $x_{3}^{2}$ | $1-x_{3}^{2}$ | 0 | $-$ | $0$
(c) | $(\sqrt{\frac{2}{3}},1,0,0)$ | All | $1$ | $0$ | $0$ | $0$ | $-1$ | $-1$
(d) | $(-\sqrt{\frac{2}{3}},-1,0,0)$ | All | $1$ | $0$ | $0$ | $0$ | $-1$ | $-1$
Table 1: Fixed points of the autonomous system.

There are four fixed points for the autonomous system of evolution equations,
which are presented in Tab.1. Fixed point (a) represents a radiation dominated
Universe, it is a saddle point since its eigenvalues are

$\mu_{(a)}=-3,2,\frac{1}{2},\frac{1}{2}.$ (39)

Fixed point (b) represents a matter dominated Universe, its eigenvalues are

$\mu_{(b)}=0,-1,\frac{-3-\sqrt{3}\sqrt{27-16x_{3}^{2}\lambda}}{4},\frac{-3+\sqrt{3}\sqrt{27-16x_{3}^{2}\lambda}}{4},$
(40)

the eigenvector corresponding to the vanishing eigenvalue reads
$(0,0,1,\frac{-x_{3}}{\sqrt{1-x_{3}^{2}}})$. Generally speaking, we need to go
to the higher order to study the stability of this fixed point, however, as
long as one can prove that one of the eigenvalues is positive, in other words,
if $\lambda<\frac{3}{2x_{3}^{2}}$ is proved to be true, such fixed point is a
saddle point. Since $x_{3}\in[0,1]$, we only need to prove that
$\lambda<\frac{3}{2}$, fortunately, this condition is satisfied, in fact, we
will show that $\lambda\ll 1$ in the next section.

Fixed point (c) and fixed point (d) are represented by two saturated de Sitter
Universes, they are always stable since their eigenvalues both read

$\displaystyle\mu_{(c)}$ $\displaystyle=$
$\displaystyle-4,-3,-\frac{3}{2},-\frac{3}{2},$ (41) $\displaystyle\mu_{(d)}$
$\displaystyle=$ $\displaystyle-4,-3,-\frac{3}{2},-\frac{3}{2}.$ (42)

One should note that the phase space is separated in two parts because of
these two symmetrical attractors, more specifically, generally speaking, the
trajectories in the phase space run toward the de Sitter attractor (c) if
$x_{2}>0$ at the beginning, otherwise they will run toward the other
attractor.

##  IV Confront the model with observations

In this section, we confront the coupled three-form dark energy model with
CMB, BAO, SN Ia observation, and OHD based on the Hubble parameter,

$H^{2}=\frac{\omega_{(r)}(1+z)^{4}+\omega_{(c)}\left(\frac{1+x_{1}^{2}}{1+x_{10}^{2}}\right)^{\frac{\lambda}{2}}(1+z)^{3}+\omega_{(b)}(1+z)^{3}}{1-x_{2}^{2}}.$
(43)

###  IV.1 CMB measurements

There are two shift parameters, $R$ and $l_{A}$, that contain much of
information of CMB power spectrum, the former is defined as

$R=\sqrt{\omega_{(m)}}(1+z_{\ast})d_{A}(z_{\ast}),$ (44)

and the latter reads

$l_{A}=\pi(1+z_{\ast})\frac{d_{A}(z_{\ast})}{r_{s}(z_{\ast})},$ (45)

where $d_{A}(z_{\ast})=\frac{1}{1+z}\int_{0}^{z_{\ast}}\frac{d\tilde{z}}{H}$
is the angular distance at decoupling Efstathiou2010Cosmic , which depends on
the dominant components after decoupling. The redshift at decoupling
$z_{\ast}$ is given by Hu1996Small

$\displaystyle z_{*}$ $\displaystyle=$ $\displaystyle
1048(1+0.00124\omega_{(b)}^{-0.738})(1+g_{1}\omega_{(m)}^{g_{2}}),$ (46)
$\displaystyle g_{1}$ $\displaystyle=$
$\displaystyle\frac{0.0783\omega_{(b)}^{-0.238}}{1+39.5\omega_{(b)}^{0.763}},$
(47) $\displaystyle g_{2}$ $\displaystyle=$
$\displaystyle\frac{0.56}{1+21.1\omega_{(b)}^{1.81}},$ (48)

where $\omega_{(m)}=\omega_{(c)}+\omega_{(b)}$. In this work, we use the
following Planck 2018 compressed likelihood chen2019distance  with these two
shift parameters to perform a likelihood analysis,

$\displaystyle\chi_{CMB}^{2}$ $\displaystyle=$ $\displaystyle
s^{T}C_{CMB}^{-1}s,$ (49) $\displaystyle s$ $\displaystyle=$
$\displaystyle(R-1.7502,l_{A}-301.471,\omega_{(b)}-0.02236),$ (50)

where $C_{ij}=D_{ij}\sigma_{i}\sigma_{j}$ is the covariance matrix,
$\sigma=(0.0046,0.09,0.00015)$ is the errors, and
$D_{CMB}=\left(\begin{array}[]{ccc}1&0.46&-0.66\\\ 0.46&1&-0.33\\\
-0.66&-0.33&1\\\ \end{array}\right)$ is the covariance.

###  IV.2 Baryon acoustic oscillations

The relative BAO distance is defined as

$d_{z}(z)=\frac{r_{s}(z_{d})}{D_{V}(z)},$ (51)

with $D_{V}(z)=[\frac{z(1+z)^{2}d_{A}(z)^{2}}{H(z)}]^{\frac{1}{3}}$. $z_{d}$
is the redshift at the drag epoch, a epoch when baryons are released from the
Compton drag of the photons. $z_{d}$ can be calculated by using
Eisenstein1997Baryonic

$\displaystyle z_{d}$ $\displaystyle=$ $\displaystyle
1291\frac{\omega_{(m)}^{0.251}}{1+0.659\omega_{(m)}^{0.828}}(1+b_{1}\omega_{(b)}^{b_{2}}),$
(52) $\displaystyle b_{1}$ $\displaystyle=$ $\displaystyle
0.313\omega_{(m)}^{-0.419}(1+0.607\omega_{(m)}^{0.674}),$ (53) $\displaystyle
b_{2}$ $\displaystyle=$ $\displaystyle 0.238\omega_{(m)}^{0.223}.$ (54)

We use four measurements from 6dFGS at $z_{eff}$= 0.106, the recent SDSS main
galaxy (MGS) at $z_{eff}$ = 0.15 ross2015clustering , and $z_{eff}$= 0.32 and
0.57 for the Baryon Oscillation Spectroscopic Survey (BOSS). Anderson2013The
.111Although the BOSS DR12 sample alam2017clustering  is available now,
fitting results won’t be too much different if we only use the old data.
Therefore, the BAO likelihood is written as

$\displaystyle\chi_{BAO}^{2}$ $\displaystyle=$ $\displaystyle
t^{T}C_{BAO}^{-1}t,$ (55) $\displaystyle t$ $\displaystyle=$
$\displaystyle(d_{z}(0.106)-0.336,d_{z}(0.15)-0.2239,d_{z}(0.32)-0.1181,d_{z}(0.57)-0.0726),$
(56)

with $C_{BAO}^{-1}=diag(4444.44,14067.88,183275.76,2004606.79)$ is the the
covariance matrix.

###  IV.3 Type Ia supernovae

For supernovae data, we employ the Joint Light-curve Analysis (JLA) sample
betoule2014improved ,222For the same reason with giving up new BAO data, we
adopt JLA sample in this paper, although the Pantheon sample
scolnic2018complete  is available now. the distance modulus is then assumed as

$\mu_{obs}=m_{B}-(M_{B}-\alpha X_{1}+\beta c),$ (57)

where $m_{B}$ and $M_{B}$ are SN Ia peak apparent magnitude and SN Ia absolute
magnitude, respectively. $\alpha$ and $\beta$ are two constants. $c$ is the
color parameter, and $X_{1}$ is the stretch factor. Therefore, the likelihood
for SN Ia is defined as

$\chi_{SN}^{2}=\Delta^{T}C_{SN}^{-1}\Delta,$ (58)

where $\Delta=\mu-\mu_{obs}$ and $C_{SN}$ is the covariance matrix.

###  IV.4 OHD

For the OHD in Table 2, the best-fit values of the model parameters can be
determined by a likelihood analysis based on the calculation of

$\chi^{2}_{H}=\sum_{i=1}^{29}\frac{(H(z_{i})-H_{obs}(z_{i}))^{2}}{\sigma_{H}^{2}(z_{i})}.$
(59)

$z$ | $H(z)$ | Ref.
---|---|---
$0.0708$ | $69.0\pm 19.68$ | Zhang et al. (2014)-zhang2014four
$0.09$ | $69.0\pm 12.0$ | Jimenez et al. (2003)-Jimenez2003
$0.12$ | $68.6\pm 26.2$ | Zhang et al. (2014)-zhang2014four
$0.17$ | $83.0\pm 8.0$ | Simon et al. (2005)-Simon2005
$0.179$ | $75.0\pm 4.0$ | Moresco et al. (2012)-Moresco2012
$0.199$ | $75.0\pm 5.0$ | Moresco et al. (2012)-Moresco2012
$0.20$ | $72.9\pm 29.6$ | Zhang et al. (2014)-zhang2014four
$0.27$ | $77.0\pm 14.0$ | Simon et al. (2005)-Simon2005
$0.28$ | $88.8\pm 36.6$ | Zhang et al. (2014)-zhang2014four
$0.352$ | $83.0\pm 14.0$ | Moresco et al. (2012)-Moresco2012
$0.3802$ | $83.0\pm 13.5$ | Moresco et al. (2016)-Moresco2016
$0.4$ | $95\pm 17.0$ | Simon et al. (2005)-Simon2005
$0.4004$ | $77.0\pm 10.2$ | Moresco et al. (2016)-Moresco2016
$0.4247$ | $87.1\pm 11.2$ | Moresco et al. (2016)-Moresco2016
$0.4497$ | $92.8\pm 12.9$ | Moresco et al. (2016)-Moresco2016
$0.4783$ | $80.9\pm 9.0$ | Moresco et al. (2016)-Moresco2016
$0.48$ | $97.0\pm 62.0$ | Stern et al. (2010)-Stern2010
$0.593$ | $104.0\pm 13.0$ | Moresco et al. (2012)-Moresco2012
$0.68$ | $92.0\pm 8.0$ | Moresco et al. (2012)-Moresco2012
$0.875$ | $125.0\pm 17.0$ | Moresco et al. (2012)-Moresco2012
$0.88$ | $90.0\pm 40.0$ | Stern et al. (2010)-Stern2010
$0.9$ | $117.0\pm 23.0$ | Simon et al. (2005)-Simon2005
$1.037$ | $154.0\pm 20.0$ | Moresco et al. (2012)-Moresco2012
$1.3$ | $168.0\pm 17.0$ | Simon et al. (2005)-Simon2005
$1.363$ | $160.0\pm 33.6$ | Moresco (2015)-Moresco2015
$1.43$ | $177.0\pm 18.0$ | Simon et al. (2005)-Simon2005
$1.53$ | $140.0\pm 14.0$ | Simon et al. (2005)-Simon2005
$1.75$ | $202.0\pm 40.0$ | Simon et al. (2005)-Simon2005
$1.965$ | $186.5\pm 50.4$ | Moresco (2015)-Moresco2015
Table 2:  The current available OHD dataset.

Finally, we have the total likelihood as

$\chi_{tot}^{2}=\chi_{CMB}^{2}+\chi_{BAO}^{2}+\chi_{SN}^{2}+\chi_{H}^{2}$ (60)

Since the probability distribution of the parameter $x_{10}$ is approximate to
a uniform distribution and its minimum approaches to minus infinity, we set
the parameter $x_{10}$ to be fixed at $-1000000$. Therefore, we choose
$\Omega_{(m)}$, $\lambda$, $h$, and $\omega_{(b)}$ as fitting parameters for
the coupled three-form dark energy model, and $\Omega_{(m)}$, $h$, and
$\omega_{(b)}$ as fitting parameters for $\Lambda$CDM. For comparison, we add
the $\Lambda$CDM model into our likelihood analysis. The fitting results are
presented in Tab.3 and Fig.1. From Tab.3, we obtain
$h=0.678_{-0.006}^{+0.007}$($1\sigma$ level) km/s/Mpc for the coupled three-
form dark energy model, and $h=0.676_{-0.005}^{+0.005}$($1\sigma$ level)
km/s/Mpc for the $\Lambda$CDM model. Corresponding, the $H_{0}$ tension
between them and $h=0.738\pm 0.011$ are reduced to $4.3\sigma$ for the coupled
three-form dark energy model, and $5.1\sigma$ for the $\Lambda$CDM model.
Therefore, within the coupled three-form dark energy model, the tension is
still strong. As shown in Fig.1, considering the extra parameters $\lambda$ in
the $\Lambda$CDM model can affect the constraints on the Hubble constant
$H_{0}$ because that $\lambda$ are positively correlated with $h$. These
fitting results are not unexpected, in fact, serval recent works on coupled
dark energy model di2020nonminimal ; di2020interacting ; cheng2020testing
produce similar results, i.e. $H_{0}$ gets a bit higher but not enough to
completely solve the $H_{0}$ tension, which is mostly alleviated by a bit
larger error bars.

We also consider the Akaike information criterion (AIC) to compare the coupled
three-form dark energy model and the $\Lambda$CDM model. The AIC is defined as
$\chi_{min}^{2}+2k$, where $k$ denotes the number of cosmological parameters.
Therefore, from Tab.3, we have AIC=$715.734$ for the coupled three-form dark
energy model and AIC=$713.873$ for the $\Lambda$CDM model. Then the
$\Lambda$CDM model is more supported than the coupled three-form dark energy
model because that, by definition, a model with a smaller value of AIC is a
more supported model.

Model | coupled three-form dark energy model |  $\Lambda$CDM
---|---|---
$\Omega_{(m)}$ | $0.323_{-0.007}^{+0.008}$ | $0.321_{-0.007}^{+0.007}$
$\lambda$ | $0.001_{-0.001}^{+0.001}$ | $-$
$h$ | $0.679_{-0.007}^{+0.008}$ | $0.676_{-0.005}^{+0.005}$
$\omega_{(b)}$ | $0.0223_{-0.0002}^{+0.0002}$ | $0.0223_{-0.0001}^{+0.0001}$
$\chi_{min}^{2}$ | 707.734 | 707.873
AIC | 715.734 | 713.873
Table 3:  Fit results of cosmological parameters in the coupled three-form
dark energy model and $\Lambda$CDM using the CMB+BAO+JLA+OHD data. Figure 1:
$1\sigma$ and $2\sigma$ confidence regions and probability densities for the
parameters in the coupled three-form dark energy model.

##  V Conclusions

In this paper, a coupled three-form dark energy model is put forward to
relieve the Hubble tension. To start with, we perform a dynamical analysis
with the coupled three-form dark energy model, and obtain three fixed points,
including a saddle point representing a radiation dominated Universe, a saddle
point representing a matter dominated Universe, and a attractor representing a
saturated de Sitter Universe. Then, we confront the $\Lambda$CDM model and the
coupled three-form dark energy model with CMB, BAO, SN Ia observation, and
obtain $h=0.676_{-0.005}^{+0.005}$($1\sigma$ level) km/s/Mpc for the
$\Lambda$CDM model, which is in strong tension with latest local $h$ value at
$5.1\sigma$ level, and $h=0.679_{-0.007}^{+0.008}$($1\sigma$ level) km/s/Mpc
for the coupled three-form dark energy model, which is in strong tension with
latest local $h$ value at $4.3\sigma$ level. Therefore, within the coupled
three-form dark energy model, the strong $H_{0}$ tension is still exist. In
addition, by using AIC, we find that the $\Lambda$CDM model is more supported
than the coupled three-form dark energy model.

## Acknowledgments

The paper is partially supported by the Natural Science Foundation of China.

## References

  * (1) N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. Banday, R. Barreiro, N. Bartolo, S. Basak, et al., arXiv preprint arXiv:1807.06209 (2018)
  * (2) A.G. Riess, S. Casertano, W. Yuan, L.M. Macri, D. Scolnic, The Astrophysical Journal 876(1), 85 (2019), arXiv:1903.07603
  * (3) K.C. Wong, S.H. Suyu, G.C. Chen, C.E. Rusu, M. Millon, D. Sluse, V. Bonvin, C.D. Fassnacht, S. Taubenberger, M.W. Auger, et al., Monthly Notices of the Royal Astronomical Society (2019), arXiv:1907.04869
  * (4) D.N. Spergel, R. Flauger, R. Hlo$\check{z}$ek, Physical Review D 91(2), 023518 (2015), arXiv:1312.3313
  * (5) G.E. Addison, Y. Huang, D.J. Watts, C.L. Bennett, M. Halpern, G. Hinshaw, J.L. Weiland, The Astrophysical Journal 818(2), 132 (2016), arXiv:1511.00055
  * (6) N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. Banday, R.B. Barreiro, N. Bartolo, S. Basak, et al., Astronomy & Astrophysics 607, A95 (2017), arXiv:1608.02487
  * (7) W. Cardona, M. Kunz, V. Pettorino, Journal of Cosmology and Astroparticle Physics 2017(03), 056 (2017), arXiv:1611.06088
  * (8) B. Follin, L. Knox, Monthly Notices of the Royal Astronomical Society 477(4), 4534 (2018), arXiv:1707.01175
  * (9) M. Li, X.D. Li, Y.Z. Ma, X. Zhang, Z. Zhang, Journal of Cosmology and Astroparticle Physics 2013(9), 87 (2013), arXiv:1305.5302
  * (10) Q. Huang, K. Wang, European Physical Journal C 76(9), 506 (2016), arXiv:1606.05965
  * (11) R.A. Battye, A. Moss, Physical Review Letters 112(5), 051303 (2014), arXiv:1308.5870
  * (12) J.F. Zhang, J.J. Geng, X. Zhang, Journal of Cosmology and Astroparticle Physics 2014(10), 044 (2014), arXiv:1408.0481
  * (13) J. Zhang, Y. Li, X. Zhang, Physics Letters B 740, 359 (2015), arXiv:1403.7028
  * (14) L. Feng, J. Zhang, X. Zhang, Science China-physics Mechanics and Astronomy 61(5), 050411 (2018), arXiv:1706.06913
  * (15) M. Zhao, J. Zhang, X. Zhang, Physics Letters B 779, 473 (2018), arXiv:1710.02391
  * (16) S.R. Choudhury, S. Choubey, European Physical Journal C 79(7), 557 (2019), arXiv:1807.10294
  * (17) S. Vagnozzi, Physical Review D 102(2), 023518 (2020), arXiv:1907.07569
  * (18) P. Agrawal, F.Y. Cyr-Racine, D. Pinner, L. Randall, arXiv preprint arXiv:1904.01016 (2019)
  * (19) V. Poulin, T.L. Smith, T. Karwal, M. Kamionkowski, Physical Review Letters 122(22), 221301 (2019), arXiv:1811.04083
  * (20) X. Li, A. Shafieloo, The Astrophysical Journal Letters 883(1), L3 (2019), arXiv:1906.08275
  * (21) D. Valentino, Eleonora, Linder, V. Eric, Melchiorri, Alessandro, Physical Review D 97(4), 043528 (2018), arXiv:1710.02153
  * (22) E. Di Valentino, A. Melchiorri, O. Mena, Physical Review D 96(4), 043503 (2017), arXiv:1704.08342
  * (23) W. Yang, A. Mukherjee, E. Di Valentino, S. Pan, Physical Review D 98(12), 123527 (2018), arXiv:1809.06883
  * (24) Y. Yao, Y. Yan, X. Meng, European Physical Journal C 78(2), 153 (2018), arXiv:1704.05772
  * (25) T.S. Koivisto, N.J. Nunes, Physical Review D 88(12), 123512 (2013), arXiv:1212.2541
  * (26) J.R. Ray, Journal of Mathematical Physics 13(10), 1451 (1972)
  * (27) T.S. Koivisto, N.J. Nunes, Physical Review D 80(10), 103509 (2009), arXiv:0908.0920
  * (28) G. Efstathiou, J.R. Bond, Monthly Notices of the Royal Astronomical Society 304(1), 75 (2010), arXiv:astro-ph/9807103
  * (29) W. Hu, N. Sugiyama, Physics 471(2), 542 (1996), arXiv:astro-ph/9510117
  * (30) L. Chen, Q.G. Huang, K. Wang, Journal of Cosmology and Astroparticle Physics 2019(02), 028 (2019), arXiv:1808.05724
  * (31) D.J. Eisenstein, W. Hu, Astrophysical Journal 496(2), 605 (1997), arXiv:astro-ph/9709112
  * (32) A.J. Ross, L. Samushia, C. Howlett, W.J. Percival, A. Burden, M. Manera, Monthly Notices of the Royal Astronomical Society 449(1), 835 (2015), arXiv:1409.3242
  * (33) L. Anderson, r. Aubourg, S. Bailey, F. Beutler, V. Bhardwaj, M. Blanton, A.S. Bolton, J. Brinkmann, J.R. Brownstein, A. Burden, Monthly Notices of the Royal Astronomical Society 441(1), 24 (2013), arXiv:1312.4877
  * (34) S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J.A. Blazek, A.S. Bolton, J.R. Brownstein, A. Burden, C.H. Chuang, et al., Monthly Notices of the Royal Astronomical Society 470(3), 2617 (2017), arXiv:1312.4877
  * (35) M.e.a. Betoule, R. Kessler, J. Guy, J. Mosher, D. Hardin, R. Biswas, P. Astier, P. El-Hage, M. Konig, S. Kuhlmann, et al., Astronomy & Astrophysics 568, A22 (2014), arXiv:1401.4064
  * (36) D.M. Scolnic, D. Jones, A. Rest, Y. Pan, R. Chornock, R. Foley, M. Huber, R. Kessler, G. Narayan, A. Riess, et al., The Astrophysical Journal 859(2), 101 (2018), arXiv:1710.00845
  * (37) C. Zhang, H. Zhang, S. Yuan, S. Liu, T.J. Zhang, Y.C. Sun, Research in Astronomy and Astrophysics 14(10), 1221 (2014), arXiv:1207.4541
  * (38) R. Jimenez, L. Verde, T. Treu, D. Stern, The Astrophysical Journal 593(2), 622 (2003), arXiv:astro-ph/0302560
  * (39) J. Simon, L. Verde, R. Jimenez, Physical Review D 71(12), 123001 (2005), arXiv:astro-ph/0412269
  * (40) M. Moresco, L. Verde, L. Pozzetti, R. Jimenez, A. Cimatti, Journal of Cosmology and Astroparticle Physics 2012(07), 053 (2012), arXiv:1201.6658
  * (41) M. Moresco, L. Pozzetti, A. Cimatti, R. Jimenez, C. Maraston, L. Verde, D. Thomas, A. Citro, R. Tojeiro, D. Wilkinson, Journal of Cosmology and Astroparticle Physics 2016(05), 014 (2016), arXiv:1601.01701
  * (42) D. Stern, R. Jimenez, L. Verde, M. Kamionkowski, S.A. Stanford, Journal of Cosmology and Astroparticle Physics 2010(02), 008 (2010), arXiv:0907.3149
  * (43) M. Moresco, Monthly Notices of the Royal Astronomical Society: Letters 450(1), L16 (2015), arXiv:1503.01116
  * (44) E. Di Valentino, A. Melchiorri, O. Mena, S. Vagnozzi, Physical Review D 101(6), 063502 (2020), arXiv:1910.09853
  * (45) E. Di Valentino, A. Melchiorri, O. Mena, S. Vagnozzi, Physics of the Dark Universe 30, 100666 (2020), arXiv:1908.04281
  * (46) G. Cheng, Y.Z. Ma, F. Wu, J. Zhang, X. Chen, Physical Review D 102(4), 043517 (2020), arXiv:1911.04520


# On a Frequency-Stabilized Single Current Inverse Source Formulation

Paolo Ricci  Adrien Merlini  Francesco P. Andriulli

###### Abstract

Several strategies are available for solving the inverse source problem in
electromagnetics. Among them, many have been focusing in retrieving Love
currents by solving, after regularization, for Love’s electric and magnetic
currents. In this work we present a dual-element discretization, analysis, and
stabilization of an inverse source formulation providing Love data by solving
for only one current. This results in substantial savings and allows for an
effective quasi-Helmholtz projector stabilization of the resulting operator.
Theoretical considerations are complemented by numerical tests showing
effectiveness and efficiency of the newly proposed method.

On a Frequency-Stabilized Single Current Inverse Source Formulation

${}^{(\ref{})}$

(1) Politecnico di Torino, Turin, Italy

(2) IMT Atlantique, Brest, France

##  I Introduction

Inverse source strategies based on external measurements are widely used for
several applications including antenna diagnostics and near-field (NF) to far-
field (FF) transformations. These methods are challenging because of the ill-
posedness of the radiation operator. The different approaches present in
literature can be divided into two categories: methods that consider both
electric and magnetic currents and deal with the null-space introduced by
enforcing an additional constraint to the solution—often the Love (or zero-
field inside) condition [1]—and strategies that partially reduce the space of
solutions by restricting the unknowns to a single type of current [1, 2]. In
both cases, the Love conditions can be enforced. While the Love constraints
have no impact on the reconstruction of the FF, their enforcement can be
useful for diagnostic purposes and can be performed either explicitly or
implicitly by leveraging an approximation of the Steklov-Poincaré operator
[2].

In this work we propose a scheme falling in the second category that targets
the reconstruction of only one of the two Love currents. Differently from
previous works, however, we propose a single source approach that does not
rely on any approximate impedance condition, but directly targets a stable
discretization of the Steklov-Poincare operator that leverages on dual
elements. This has two advantages: in addition to providing Love currents
without scattering approximations, it enables the stabilization of the
formulation at lower frequencies. A second contribution of this work is, in
fact, a low-frequency stabilization of the Steklov-Poincaré operator based on
quasi-Helmholtz projectors. The effectiveness of the new schemes is validated
both by theoretical considerations and by numerical results that show the
relevance of the methods in real case scenarios.

##  II Notation and Background

Let $\Gamma$ be a closed, simply connected surface surrounding an electric
source. By means of the equivalence theorem, equivalent electric and magnetic
surface currents densities $\boldsymbol{J}$ and $\boldsymbol{M}$ can be
obtained on $\Gamma$ that radiate the same external field as the original
source. We define the operators

$\displaystyle{\mathcal{T}}(\boldsymbol{J})(\boldsymbol{r})$
$\displaystyle=k\,\mathcal{T}_{s}(\boldsymbol{J})(\boldsymbol{r})+k^{-1}\,\mathcal{T}_{h}(\boldsymbol{J})(\boldsymbol{r})\,,$
(1) $\displaystyle{\mathcal{K}}(\boldsymbol{J})(\boldsymbol{r})$
$\displaystyle=-\boldsymbol{\hat{n}_{{r}}}\times
p.v.\int_{\Gamma}\nabla\times\frac{e^{ik\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|}}{4\pi\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|}\boldsymbol{J}(\boldsymbol{r}^{\prime})\,\mathrm{d}s^{\prime}$
(2)

with
$\mathcal{T}_{s}(\boldsymbol{J})(\boldsymbol{r})=i\boldsymbol{\hat{n}_{{r}}}\times\int_{\Gamma}\frac{e^{ik\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|}}{4\pi\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|}{\boldsymbol{J}(\boldsymbol{r}^{\prime})}\,\mathrm{d}\boldsymbol{r}^{\prime}$,
$\mathcal{T}_{h}(\boldsymbol{J})(\boldsymbol{r})=i\boldsymbol{\hat{n}_{{r}}}\times\nabla\int_{\Gamma}\frac{e^{ik\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|}}{4\pi\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|}{\nabla_{s}\cdot\boldsymbol{J}(\boldsymbol{r}^{\prime})}\,\mathrm{d}\boldsymbol{r}^{\prime}$,
and where $k$ is the wavenumber and $\boldsymbol{\hat{n}_{{r}}}$ is outside
pointing normal vector at $\boldsymbol{r}\in\Gamma$. Let ${\mathcal{P}}^{+}$
and ${\mathcal{P}}^{-}$ be respectively the external and the internal Calderón
projectors that satisfy the property
${\mathcal{P}}^{-}{\mathcal{P}}^{+}=\boldsymbol{0}$, from which we obtain that

${\mathcal{P}}^{-}\begin{bmatrix}-\boldsymbol{M}^{e}\\\
\eta_{0}\boldsymbol{J}^{e}\end{bmatrix}=\begin{bmatrix}\frac{{\mathcal{I}}}{2}+{\mathcal{K}}&-{\mathcal{T}}\\\
{\mathcal{T}}&\frac{{\mathcal{I}}}{2}+{\mathcal{K}}\end{bmatrix}\begin{bmatrix}-\boldsymbol{M}^{e}\\\
\eta_{0}\boldsymbol{J}^{e}\end{bmatrix}=\boldsymbol{0}$ (3)

where $\boldsymbol{M}^{e}$ and $\boldsymbol{J}^{e}$ are Love magnetic and
electric currents on $\Gamma$ (radiating zero-fields inside and the original
field outside). Given a second measurement surface $\Gamma_{m}$ external to
$\Gamma$, we denote by

${\mathcal{R}}=\begin{bmatrix}-{\mathcal{K}}_{\boldsymbol{r}}&{\mathcal{T}}_{\boldsymbol{r}}\\\
-{\mathcal{T}}_{\boldsymbol{r}}&-{\mathcal{K}}_{\boldsymbol{r}}\end{bmatrix}$
(4)

the tangential projection of the radiation operator that defines the relation
between the currents on $\Gamma$ and the tangential fields on the measurement
surface $\Gamma_{m}$, with ${\mathcal{T}}_{\boldsymbol{r}}={\mathcal{T}}$. The
operator ${\mathcal{K}}_{\boldsymbol{r}}$ is similar to ${\mathcal{K}}$
without the principal value. The corresponding linear system is

${\mathcal{R}}\begin{bmatrix}-\boldsymbol{M}\\\
\eta_{0}\boldsymbol{J}\end{bmatrix}=\begin{bmatrix}-\boldsymbol{E}\\\
\eta_{0}\boldsymbol{H}\end{bmatrix}\,.$ (5)

##  III A Stable Discretization of the Steklov-Poincaré Operator

In a general inverse source setting, current sources can be recovered by
pseudo-inverting (or iteratively pseudo-inverting) the discretized counterpart
of ${\mathcal{R}}$ in (5) and to which additional equations are added to
enforce the Love condition [1]. In this work, instead, we consider a single
row of (4) (e.g. the first row, related to the electric field) and we formally
solve for one of the currents which yields the equation

$\left(-{\mathcal{K}}_{\boldsymbol{r}}-{\mathcal{T}}_{\boldsymbol{r}}\left(\frac{{\mathcal{I}}}{2}+{\mathcal{K}}\right)^{-1}{\mathcal{T}}\right)\left(\boldsymbol{M}\right)=\boldsymbol{E}\,,$
(6)

where $-\left(\frac{{\mathcal{I}}}{2}+{\mathcal{K}}\right)^{-1}{\mathcal{T}}$
is the Steklov-Poincaré operator. Equation (6) can, upon discretization, be
pseudo-inverted using a Moore-Penrose pseudo-inverse to obtain the magnetic
Love current. In practice, however, the stable discretization of (6) cannot be
achieved by only using the well known Rao-Wilton-Glisson (RWG) functions. One
possibility is to rely on an approximation of the Steklov-Poincaré operator as
is done in [2]. Here we will follow a different approach based on dual
elements which, in addition to providing a consistent discretization of the
operator, will also allow for its stabilization with quasi-Helmholtz
projectors as will be discussed in the next section. After approximating the
geometries of $\Gamma$ and of $\Gamma_{m}$ with meshes of triangular elements
with $N_{e}$ and $N_{m}$ edges respectively, the discretization we propose is

$\left(-\mathbf{{K}}_{\boldsymbol{r}}-\mathbb{{T}}_{\boldsymbol{r}}\left(\mathbb{{G}}/2+\mathbb{{K}}\right)^{-1}\mathbf{{T}}\right){\mathbf{m}}={\mathbf{e}}$
(7)

where
$\boldsymbol{E}(\boldsymbol{r})\approx\sum_{m=1}^{N_{m}}\left[{\mathbf{e}}\right]_{m}\boldsymbol{\hat{n}_{{r}}}\times\boldsymbol{g}_{m}(\boldsymbol{r})$,
$\boldsymbol{M}(\boldsymbol{r})\approx\sum_{n=1}^{N_{e}}\left[\boldsymbol{m}\right]_{n}\boldsymbol{f}_{n}(\boldsymbol{r})$,
$\boldsymbol{f}(\boldsymbol{r})$ are the RWG basis functions (defined without
edge normalization), $\boldsymbol{g}(\boldsymbol{r})$ are the Buffa-
Christiansen basis functions (see [3] for their definition), and where
$[\mathbf{{K}}_{\boldsymbol{r}}]_{ij}=\left<\boldsymbol{\hat{n}_{{r}}}\times\boldsymbol{g}_{i},{\mathcal{K}}_{\boldsymbol{r}}\boldsymbol{f}_{j}\right>_{\Gamma_{m}}$,
$[\mathbb{{T}}_{\boldsymbol{r}}]_{ij}=\left<\boldsymbol{\hat{n}_{{r}}}\times\boldsymbol{g}_{i},{\mathcal{T}}\boldsymbol{g}_{j}\right>_{\Gamma_{m}}$,
$[\mathbb{{G}}]_{ij}=\left<\boldsymbol{\hat{n}_{{r}}}\times\boldsymbol{f}_{i},\boldsymbol{g}_{j}\right>_{\Gamma}$,
$[\mathbb{{K}}]_{ij}=\left<\boldsymbol{\hat{n}_{{r}}}\times\boldsymbol{f}_{i},{\mathcal{K}}\boldsymbol{g}_{j}\right>_{\Gamma}$,
and
$[\mathbf{{T}}]_{ij}=\left<\boldsymbol{\hat{n}_{{r}}}\times\boldsymbol{f}_{i},{\mathcal{T}}\boldsymbol{f}_{j}\right>_{\Gamma}$,
with
$\left<\boldsymbol{a},\boldsymbol{b}\right>_{\Gamma}=\int_{\Gamma}\boldsymbol{a}(\boldsymbol{r})\cdot\boldsymbol{b}(\boldsymbol{r})\,\mathrm{d}s$.
After solving (7) and recovering ${\mathbf{m}}$, the other Love current can be
obtained by applying the Steklov-Poincaré operator once more.

##  IV Low-Frequency Stabilization

The operator in (7) (as well as the original system (5)) is affected by low
frequency breakdown, i.e. its condition number increases as the frequency
decreases, leading to the insurgence of numerical errors. The stabilization we
propose is based on Quasi-Helmoltz projectors that are defined as
$\mathbf{{P}}^{\Sigma}=\mathbf{{\Sigma}}(\mathbf{{\Sigma}}^{\text{T}}\mathbf{{\Sigma}})^{+}\mathbf{{\Sigma}}^{\text{T}}$,
$\mathbf{{P}}^{\Lambda H}=\mathbf{{I}}-\mathbf{{P}}^{\Sigma}$,
$\mathbb{{P}}^{\Lambda}=\mathbf{{\Lambda}}(\mathbf{{\Lambda}}^{\text{T}}\mathbf{{\Lambda}})^{+}\mathbf{{\Lambda}}^{\text{T}}$,
and $\mathbb{{P}}^{\Sigma H}=\mathbf{{I}}-\mathbb{{P}}^{\Lambda}$, where
$\mathbf{{\Lambda}}$ and $\mathbf{{\Sigma}}$ are respectively the loop- and
the star-decomposition matrices defined in [4], and + denotes the Moore-
Penrose pseudo-inverse. The stabilized equation reads

$\mathbb{{M}}\left(-\mathbf{{K}}_{\boldsymbol{r}}-(\mathbb{{T}}_{\boldsymbol{r}}\left(\mathbb{{G}}/2+\mathbb{{K}}\right)^{-1}\mathbf{{T}}\right)\mathbf{{M}}{\mathbf{x}}=\mathbb{{M}}{\mathbf{e}}$
(8)

where $\mathbf{{M}}{\mathbf{x}}={\mathbf{m}}$,
$\mathbf{{M}}=1/\sqrt{k}\mathbf{{P}}^{\Lambda
H}+i\sqrt{k}\mathbf{{P}}^{\Sigma}$ and $\mathbb{{M}}$ is the same as
$\mathbf{{M}}$ with the roles of $\mathbf{{\Lambda}}$ and $\mathbf{{\Sigma}}$
exchanged. We have to omit the technical details for the sake of conciseness,
but we can rigorously prove the stabilization by leveraging the cancellations
of ${\mathcal{T}}_{h}$ on the solenoidal subspaces, the mutual orthogonality
of loop and star coefficients, and finally on our proof of the property
$\lim_{k\to
0}\|\mathbb{{P}}^{\Lambda}\left(\mathbb{{G}}/2+\mathbb{{K}}\right)^{-1}\mathbf{{P}}^{\Sigma}\|=0$.

##  V Numerical Results

The effectiveness of the approach we proposed has been verified by imaging an
electric dipole oscillating at $3.16\text{\,}\mathrm{GHz}$ located in the
center of a spherical equivalent surface $\Gamma$ with a radius of
$4\text{\,}\mathrm{cm}$. First, the equivalent currents have been
reconstructed from the fields scattered by the dipole at a distance of one
wavelength from $\Gamma$; these currents have then been radiated on several
concentric spherical surfaces with increasing radii. The error of the
reconstructed fields with respect to the ones obtained from the dipole is
depicted in Fig. 1 and clearly show that the Steklov-Poincaré-based approach
achieves similar results in the reconstruction of the fields with respect to
the popular standard formulation—based on the solution of (5) with Love
constraints—despite being involving a significantly smaller linear system.
Furthermore we verified that the condition number of the stabilized operator
evaluated at the pseudo-inversion threshold does not increase when the
frequency decreases (Fig. 2), showing that the new method is free from low
frequency conditioning breakdown.

Figure 1: Relative error of the reconstruction of the fields over the distance
from the equivalent surface in terms of wavelength. Figure 2: Condition number
of the stabilized operator in function of the frequency.

## ACKNOWLEDGEMENT

This work was supported in part by the European Research Council (ERC) through
the European Union’s Horizon 2020 Research and Innovation Programme under
Grant 724846 (Project 321) and in part by the H2020-MSCA-ITN-EID project
COMPETE GA No 955476.

## References

  * [1] J. L. A. Quijano and G. Vecchi, “Field and source equivalence in source reconstruction on 3d surfaces,” _Progress In Electromagnetics Research_ , vol. 103, pp. 67–100, 2010.
  * [2] J. Kornprobst, J. Knapp, R. A. Mauermayer, O. Neitz, A. Paulus, and T. F. Eibert, “Accuracy and conditioning of surface-source based near- field to far-field transformations,” _IEEE Transactions on Antennas and Propagation_ , vol. 69, no. 8, pp. 4894–4908, 2021.
  * [3] F. P. Andriulli, K. Cools, H. Bagci, F. Olyslager, A. Buffa, S. Christiansen, and E. Michielssen, “A multiplicative calderon preconditioner for the electric field integral equation,” _IEEE Transactions on Antennas and Propagation_ , vol. 56, no. 8, pp. 2398–2412, 2008.
  * [4] F. P. Andriulli, K. Cools, I. Bogaert, and E. Michielssen, “On a well- conditioned electric field integral operator for multiply connected geome- tries,” _IEEE transactions on antennas and propagation_ , vol. 61, no. 4, pp. 2077–2087, 2012.