Observability of the Higgs boson decay to a photon and a dark photon

We begin by introducing the relevant fields and Lagrangian. Consider a vector-like fermion $\psi_{1}$ that transforms under a representation of $SU(2)_{L}$ of dimension $p=n\pm 1$, has a weak hypercharge of $Y^{p}=Y^{n}+1/2$ and a charge $Q^{\prime}$ under $U(1)^{\prime}$. Consider another vector-like fermion $\psi_{2}$ that transforms under a representation of $SU(2)_{L}$ of dimension $n$, has a weak hypercharge of $Y^{n}$, and a charge $Q^{\prime}$ under $U(1)^{\prime}$. The Lagrangian that determines the masses of the fermions is

	$\displaystyle\mathcal{L}_{m}=$	$\displaystyle-\left[\sum_{a,b,c}\hat{d}^{pn}_{abc}\overline{\psi}_{1}^{a}(A_{L}P_{L}+A_{R}P_{R})\psi_{2}^{b}H^{c}+\text{H.c.}\right]$		(2)
		$\displaystyle-\mu_{1}\overline{\psi}_{1}\psi_{1}-\mu_{2}\overline{\psi}_{2}\psi_{2}.$

where $a$, $b$, and $c$ are $SU(2)_{L}$ indices and are summed from 1 (corresponding to the highest weight state) to the size of the corresponding multiplet. We will write explicitly $SU(2)_{L}$ indices and sums when they are non-trivial. The $SU(2)_{L}$ tensor $\hat{d}^{pn}_{abc}$ is uniquely fixed by gauge invariance and given by the Clebsch-Gordan coefficient

$$\hat{d}^{pn}_{abc}=C^{JM}_{j_{1}m_{1}j_{2}m_{2}}=\langle j_{1}j_{2}m_{1}m_{2}|JM\rangle,$$

(3)

where

	$\displaystyle J$	$\displaystyle=\frac{p-1}{2},$	$\displaystyle j_{1}$	$\displaystyle=\frac{n-1}{2},$	$\displaystyle j_{2}$	$\displaystyle=\frac{1}{2},$		(4)
	$\displaystyle M$	$\displaystyle=\frac{p+1-2a}{2},$	$\displaystyle m_{1}$	$\displaystyle=\frac{n+1-2b}{2},$	$\displaystyle m_{2}$	$\displaystyle=\frac{3-2c}{2}.$

We will assume throughout this paper that, as is most commonly the case, the phase convention of the Clebsch-Gordan coefficients is chosen such that they are always real. We will further assume that the Clebsch-Gordan coefficients that correspond to non-physical spin combinations are zero. These two assumptions will lighten the notation. The parameters $A_{L}$ and $A_{R}$ can be complex, and it is easy to verify that there is a complex phase that cannot generally be reabsorbed via field redefinition. Such a physical complex phase will lead to CP-violating interactions. Once the Higgs field obtains a vacuum expectation value (VEV), the Lagrangian $\mathcal{L}_{m}$ will contain the mass terms

	$\displaystyle\mathcal{L}_{m}\supset$	$\displaystyle-\Biggl{[}\sum_{a,b}\frac{A_{L}v}{\sqrt{2}}\hat{d}^{pn}_{ab2}\overline{\psi}^{a}_{1}P_{L}\psi_{2}^{b}+\sum_{a,b}\frac{A_{R}^{\ast}v}{\sqrt{2}}\hat{d}^{pn}_{ba2}\overline{\psi}_{2}^{a}P_{L}\psi_{1}^{b}$		(5)
		$\displaystyle\hskip 14.22636pt+\mu_{1}\overline{\psi}_{1}P_{L}\psi_{1}+\mu_{2}\overline{\psi}_{2}P_{L}\psi_{2}+\text{H.c.}\Biggr{]},$

where $v\approx 246$ GeV is the VEV of the Higgs field. Introduce the convenient notation

$$\hat{\psi}=\begin{pmatrix}\psi_{1}\\ \psi_{2}\end{pmatrix}$$

(6)

and

$$\quad d^{pn}_{ab}=\begin{cases}\hat{d}^{pn}_{a(b-p)2},&\text{if $a\in[1,p]$ and $b\in[p+1,p+n]$},\\ 0,&\text{otherwise.}\end{cases}\\$$

(7)

The mass Lagrangian can be written more succinctly as

$$\mathcal{L}_{m}\supset-\sum_{a,b}M_{ab}\overline{\hat{\psi}}^{a}P_{L}\hat{\psi}^{b}+\text{H.c.},$$

(8)

where the mass matrix is

$$M=\begin{pmatrix}\mu_{1}\mathbbm{1}_{p\times p}&0_{p\times n}\\ 0_{n\times p}&\mu_{2}\mathbbm{1}_{n\times n}\end{pmatrix}+\frac{A_{L}v}{\sqrt{2}}d^{pn}+\frac{A_{R}^{\ast}v}{\sqrt{2}}d^{pnT}.$$

(9)

The mass matrix can then be diagonalized by introducing the fields $\tilde{\psi}^{a}$ via

$$P_{L}\hat{\psi}=R_{L}P_{L}\tilde{\psi},\quad P_{R}\hat{\psi}=R_{R}P_{R}\tilde{\psi},$$

(10)

where $R_{L}$ and $R_{R}$ are unitary matrices that diagonalize $M^{\dagger}M$ and $MM^{\dagger}$, respectively. Of course, these matrices only mix particles of identical electric charges and all entries that correspond to mixing of particles of different charges are zero. The fields $\tilde{\psi}^{a}$ are the mass eigenstates of mass $m_{a}$, and there are $p+n$ of them.

Finally, the interactions of the Higgs boson with the mass eigenstates $\tilde{\psi}^{a}$ are controlled by the terms

$$\mathcal{L}_{m}\supset-\sum_{a,b}\Omega_{ab}h\overline{\tilde{\psi}}^{a}P_{L}\tilde{\psi}^{b}+\text{H.c.},$$

(11)

where $\Omega$ is generally neither Hermitian nor diagonal and given by

$$\quad\Omega=\frac{A_{L}}{\sqrt{2}}R_{R}^{\dagger}d^{pn}R_{L}+\frac{A_{R}^{\ast}}{\sqrt{2}}R_{R}^{\dagger}d^{pnT}R_{L}.$$

(12)

We now discuss all relevant gauge interactions of the mass eigenstates $\tilde{\psi}^{a}$. The interactions of the $A/A^{\prime}$ with $\tilde{\psi}^{a}$ are controlled by

$$\mathcal{L}_{g}\supset-eA_{\mu}\overline{\tilde{\psi}}\gamma^{\mu}\tilde{Q}\tilde{\psi}-Q^{\prime}e^{\prime}A_{\mu}^{\prime}\overline{\tilde{\psi}}\gamma^{\mu}\tilde{\psi},$$

(13)

where $e^{\prime}$ is the gauge coupling constant of $U(1)^{\prime}$ and $\tilde{Q}$ is the diagonal charge matrix

$$\tilde{Q}=R_{L}^{\dagger}\hat{Q}R_{L}=R_{R}^{\dagger}\hat{Q}R_{R},$$

(14)

as QED is a vector-like interaction, where

$$\hat{Q}=\begin{pmatrix}Y^{p}+T_{3}^{p}&0_{p\times n}\\ 0_{n\times p}&Y^{n}+T_{3}^{n}\end{pmatrix},$$

(15)

with $(T_{3}^{p})_{ab}=(p+1-2a)\delta_{ab}/2$ and similarly for $T_{3}^{n}$. In practice, $\tilde{Q}$ is identical to $\hat{Q}$ except for a potential reordering of the diagonal elements.

The interactions between the $Z$ boson and $\tilde{\psi}^{a}$ is controlled by the terms

$$\mathcal{L}_{g}\supset-\sqrt{g^{2}+{g^{\prime}}^{2}}Z_{\mu}\overline{\tilde{\psi}}\gamma^{\mu}(B_{L}P_{L}+B_{R}P_{R})\tilde{\psi},$$

(16)

where $B_{L}$ and $B_{R}$ are, in general, non-diagonal but Hermitian and given by

	$\displaystyle B_{L}$	$\displaystyle=R_{L}^{\dagger}\begin{pmatrix}-s_{W}^{2}Y^{p}+c_{W}^{2}T_{3}^{p}&0_{p\times n}\\ 0_{n\times p}&-s_{W}^{2}Y^{n}+c_{W}^{2}T_{3}^{n}\end{pmatrix}R_{L},$		(17)
	$\displaystyle B_{R}$	$\displaystyle=R_{R}^{\dagger}\begin{pmatrix}-s_{W}^{2}Y^{p}+c_{W}^{2}T_{3}^{p}&0_{p\times n}\\ 0_{n\times p}&-s_{W}^{2}Y^{n}+c_{W}^{2}T_{3}^{n}\end{pmatrix}R_{R},$

with $g$ ($g^{\prime}$) denoting the $SU(2)_{L}$ ($U(1)_{Y}$) gauge coupling constant and $s_{W}$ ($c_{W}$) the sine (cosine) of the weak angle.

The interactions between the $W$ boson and $\tilde{\psi}^{a}$ is controlled by the terms

$$\mathcal{L}_{g}\supset-\frac{g}{\sqrt{2}}\overline{\tilde{\psi}}\gamma^{\mu}\left(\hat{A}_{L}P_{L}W^{+}_{\mu}+\hat{A}_{R}P_{R}W^{+}_{\mu}\right)\tilde{\psi}+\text{H.c.},$$

(18)

where

$$\hat{A}_{L}=R_{L}^{\dagger}\begin{pmatrix}T_{+}^{p}&0_{p\times n}\\ 0_{n\times p}&T_{+}^{n}\end{pmatrix}R_{L},\;\hat{A}_{R}=R_{R}^{\dagger}\begin{pmatrix}T_{+}^{p}&0_{p\times n}\\ 0_{n\times p}&T_{+}^{n}\end{pmatrix}R_{R},$$

(19)

with $(T_{+}^{p})_{ab}=\sqrt{a(p-a)}\delta_{a,b-1}$ and similarly for $T_{+}^{n}$.

The interactions of the mediators $\tilde{\psi}^{a}$ lead to contributions to the amplitude of the Higgs decay to $AA^{\prime}$, but also to those of the experimentally constrained decays to $AA$ and $A^{\prime}A^{\prime}$. The relevant diagrams are shown in Fig. 1(a). Irrespective of the mediators, gauge invariance forces the amplitudes to take the forms

	$\displaystyle M^{h\to AA}=$	$\displaystyle S^{h\to AA}\left(p_{1}\cdot p_{2}g_{\mu\nu}-p_{1\mu}p_{2\nu}\right)\epsilon^{\nu}_{p_{1}}\epsilon^{\mu}_{p_{2}}$		(20)
		$\displaystyle+i\tilde{S}^{h\to AA}\epsilon_{\mu\nu\alpha\beta}p_{1}^{\alpha}p_{2}^{\beta}\epsilon^{\nu}_{p_{1}}\epsilon^{\mu}_{p_{2}},$
	$\displaystyle M^{h\to AA^{\prime}}=$	$\displaystyle S^{h\to AA^{\prime}}\left(p_{1}\cdot p_{2}g_{\mu\nu}-p_{1\mu}p_{2\nu}\right)\epsilon^{\nu}_{p_{1}}\epsilon^{\mu}_{p_{2}}$
		$\displaystyle+i\tilde{S}^{h\to AA^{\prime}}\epsilon_{\mu\nu\alpha\beta}p_{1}^{\alpha}p_{2}^{\beta}\epsilon^{\nu}_{p_{1}}\epsilon^{\mu}_{p_{2}},$
	$\displaystyle M^{h\to A^{\prime}A^{\prime}}=$	$\displaystyle S^{h\to A^{\prime}A^{\prime}}\left(p_{1}\cdot p_{2}g_{\mu\nu}-p_{1\mu}p_{2\nu}\right)\epsilon^{\nu}_{p_{1}}\epsilon^{\mu}_{p_{2}}$
		$\displaystyle+i\tilde{S}^{h\to A^{\prime}A^{\prime}}\epsilon_{\mu\nu\alpha\beta}p_{1}^{\alpha}p_{2}^{\beta}\epsilon^{\nu}_{p_{1}}\epsilon^{\mu}_{p_{2}},$

where $p_{1}$ and $p_{2}$ are the momenta of the two gauge bosons. The $S$ coefficients are CP conserving, while the $\tilde{S}$ coefficients are CP violating. For the fermion mediators, the coefficients are given at one loop by

$$\begin{split}S^{h\to AA}&=e^{2}\sum_{a}\text{Re}(\Omega_{aa})\tilde{Q}_{aa}^{2}S_{a}+S^{h\to AA}_{\text{SM}},\\ S^{h\to AA^{\prime}}&=ee^{\prime}\sum_{a}\text{Re}(\Omega_{aa})\tilde{Q}_{aa}Q^{\prime}S_{a},\\ S^{h\to A^{\prime}A^{\prime}}&={e^{\prime}}^{2}\sum_{a}\text{Re}(\Omega_{aa}){Q^{\prime}}^{2}S_{a},\\ \tilde{S}^{h\to AA}&=e^{2}\sum_{a}\text{Im}(\Omega_{aa})\tilde{Q}_{aa}^{2}\tilde{S}_{a}+\tilde{S}^{h\to AA}_{\text{SM}},\\ \tilde{S}^{h\to AA^{\prime}}&=ee^{\prime}\sum_{a}\text{Im}(\Omega_{aa})\tilde{Q}_{aa}Q^{\prime}\tilde{S}_{a},\\ \tilde{S}^{h\to A^{\prime}A^{\prime}}&={e^{\prime}}^{2}\sum_{a}\text{Im}(\Omega_{aa}){Q^{\prime}}^{2}\tilde{S}_{a},\end{split}$$

(21)

where $S^{h\to AA}_{\text{SM}}\approx 3.3\times 10^{-5}$ GeV^-1 and $\tilde{S}^{h\to AA}_{\text{SM}}\approx 0$ GeV^-1 are the SM contributions to their respective coefficients and

	$\displaystyle S_{a}$	$\displaystyle=\frac{-m_{a}}{2\pi^{2}m_{h}^{2}}\left[2+(4m_{a}^{2}-m_{h}^{2})C_{0}(0,0,m_{h}^{2};m_{a},m_{a},m_{a})\right],$		(22)
	$\displaystyle\tilde{S}_{a}$	$\displaystyle=-i\frac{m_{a}}{2\pi^{2}}C_{0}(0,0,m_{h}^{2};m_{a},m_{a},m_{a}),$

with $C_{0}(s_{1},s_{12},s_{2};m_{0},m_{1},m_{2})$ being the scalar three-point Passarino-Veltman function Passarino:1978jh .²²2All loop calculations in this paper are performed with the help of Package-X Patel:2015tea . The decay widths are then given by

$$\begin{gathered}\Gamma^{h\to AA}=\frac{|S^{h\to AA}|^{2}+|\tilde{S}^{h\to AA}|^{2}}{64\pi}m_{h}^{3},\\ \Gamma^{h\to AA^{\prime}}=\frac{|S^{h\to AA^{\prime}}|^{2}+|\tilde{S}^{h\to AA^{\prime}}|^{2}}{32\pi}m_{h}^{3},\\ \Gamma^{h\to A^{\prime}A^{\prime}}=\frac{|S^{h\to A^{\prime}A^{\prime}}|^{2}+|\tilde{S}^{h\to A^{\prime}A^{\prime}}|^{2}}{64\pi}m_{h}^{3},\end{gathered}$$

(23)

where a symmetry factor of $1/2$ has been included for the $AA$ and $A^{\prime}A^{\prime}$ final states. The decay of the Higgs boson to a $Z$ boson and either $A$ or $A^{\prime}$ is shown in Fig. 1(b) and has an amplitude similar in form to Eq. (21), albeit with much more complicated coefficients $S^{h\to ZA}$, $\tilde{S}^{h\to ZA}$, $S^{h\to ZA^{\prime}}$, and $\tilde{S}^{h\to ZA^{\prime}}$ due to the possibility of two different mediators running in the loop.

There are several crucial results of this section that need to be discussed. First, the presence of the Levi-Civita tensor in the amplitudes is due to the $\gamma^{5}$ Dirac matrix in the Higgs to two mediators vertex. There will be no such term for the scalar cases.

Second, the amplitudes of Eq. (20) are all highly correlated. Because of this, a large $\text{BR}(h\to AA^{\prime})$ will generally lead to either a large branching ratio of the Higgs boson to invisible particles, a large modification of the coupling of the Higgs to photons, or both. The Higgs signal strengths will therefore impose strong constraints on $\text{BR}(h\to AA^{\prime})$.

Third, the decay width into two photons will in general contain terms that come from the interference between the SM amplitudes and the mediator amplitudes. When present, these cross terms typically dominate the modifications of the decay width. The presence of cross terms could, in principle, be avoided in two ways. The first way would be for the mediators to only provide a purely imaginary contribution to $S^{h\to AA}$. Since the SM contribution to $S^{h\to AA}$ is almost purely real, there would be essentially no interference. However, this turns out to be impossible for the fermion mediators. As seen in Eq. (21), all constants appearing in the $S^{h\to AA}$ are real. In addition, the kinematic function $S_{a}$ must be purely real, as bounds from the Large Electron-Positron collider (LEP) prevent charged mediators from being sufficiently light to be able to “cut” the diagram LEP1 ; LEP2 . The second way to avoid interference terms would be for the mediators to only contribute to $\tilde{S}^{h\to AA}$. This can indeed be done for fermion mediators and the signal strength constraints can therefore mostly be evaded. However, a large $\text{BR}(h\to AA^{\prime})$ will in this case lead to a large EDM for the electron. The limits on the EDM will force the complex phase to be small and will close this loophole.

Fourth, it is possible at this point to perform a naive estimate of the upper limit on $\text{BR}(h\to AA^{\prime})$ allowed by the Higgs signal strengths. A sufficiently large $\text{BR}(h\to AA^{\prime})$ will imply large Yukawa couplings. This will generally split the masses of the different $\tilde{\psi}^{a}$ and lead to a single mediator dominating the amplitude. Assume as justified above that $\text{Im}(\Omega_{ii})=\text{Im}(S_{a})=0$. Call $\Delta\text{BR}(h\to AA)$ the deviation of $\text{BR}(h\to AA)$ from its SM value and assume that it is small. Then, the following approximate relation holds:

		$\displaystyle\text{BR}(h\to AA^{\prime})\approx$		(24)
		$\displaystyle\hskip 14.22636pt\sqrt{\text{BR}(h\to A^{\prime}A^{\prime})\text{BR}(h\to AA)}\left|\frac{\Delta\text{BR}(h\to AA)}{\text{BR}(h\to AA)}\right|.$

The branching ratio $\text{BR}(h\to AA)$ is about $0.23\%$ and can at most deviate by $\mathcal{O}(25\%)$ from this value. The branching ratio of the Higgs to invisible particles $\text{BR}(h\to A^{\prime}A^{\prime})$ is at most $\mathcal{O}(10\%)$. This means that $\text{BR}(h\to AA^{\prime})$ is at most $\mathcal{O}(0.4\%)$. We will see that this approximation holds well, though other constraints will often force $\text{BR}(h\to AA^{\prime})$ to be even smaller. Of course, Eq. (24) assumes that the amplitude is dominated by a single mediator, which may not be a good approximation in more complicated models. Deviations will be observed when multiple mediators contribute comparably to $\text{BR}(h\to AA^{\prime})$, but the limits will remain of the same order of magnitude barring large fine-tuning or elaborate model building.

The constraints associated with the Higgs signal strengths are taken into account by using the $\kappa$ formalism Heinemeyer:2013tqa . Assume a production mechanism $i$ with cross section $\sigma_{i}$ or a decay process $i$ with width $\Gamma_{i}$. The parameter $\kappa_{i}$ is defined such that

$$\kappa_{i}^{2}=\frac{\sigma_{i}}{\sigma_{i}^{\text{SM}}}\quad\text{or}\quad\kappa_{i}^{2}=\frac{\Gamma_{i}}{\Gamma_{i}^{\text{SM}}},$$

(25)

where $\sigma_{i}^{\text{SM}}$ and $\Gamma_{i}^{\text{SM}}$ are the corresponding SM quantities. The only two Higgs couplings that are affected at leading order are those associated with $AA$ and $AZ$. The corresponding $\kappa$’s are

	$\displaystyle\kappa_{AA}^{2}=\frac{|S^{h\to AA}|^{2}+|\tilde{S}^{h\to AA}|^{2}}{|S^{h\to AA}_{\text{SM}}|^{2}+|\tilde{S}^{h\to AA}_{\text{SM}}|^{2}},$		(26)
	$\displaystyle\kappa_{ZA}^{2}=\frac{|S^{h\to ZA}|^{2}+|\tilde{S}^{h\to ZA}|^{2}}{|S^{h\to ZA}_{\text{SM}}|^{2}+|\tilde{S}^{h\to ZA}_{\text{SM}}|^{2}}.$

The invisible ($A^{\prime}A^{\prime}$) and semi-invisible ($A^{\prime}A$ and $A^{\prime}Z$) decays of the Higgs boson are taken into account by properly rescaling the signal strengths. The resulting global reduction of the Higgs signal strengths renders constraints from searches for the Higgs decay to invisible particles mostly superfluous and we do not impose them ATLAS:2023tkt . Constraints are then applied by using the Higgs signal strength measurements of Ref. CMS-PAS-HIG-19-005 by CMS, which uses $137$ fb^-1 of integrated luminosity at 13 TeV center-of-mass energy, and Ref. ATLAS-CONF-2021-053 by ATLAS, which uses $139$ fb^-1 of integrated luminosity at also 13 TeV. These studies conveniently provide all the information necessary (measurements, uncertainties, and correlations) to produce our own $\chi^{2}$ fit. The two searches are assumed to be uncorrelated. As we will be interested in two-dimensional scans, a point of parameter space will be considered excluded at 95% CL if its $\chi^{2}$ satisfies $\chi^{2}-\chi^{2}_{\text{min}}>5.99$, where $\chi^{2}_{\text{min}}$ is the best fit of the model.³³3Technically speaking, it is possible for the mediator contributions to $S^{h\to AA}$ to be about $-2S^{h\to AA}_{\text{SM}}$. This would mostly avoid the bounds from the Higgs signal strengths and could in principle lead to a larger $\text{BR}(h\to AA^{\prime})$. It however requires very exotic circumstances and a very large amount of fine-tuning. As such, we will not consider such cases.

As explained in Sec. III.3, most constraints from the Higgs signal strengths can be evaded by having almost purely imaginary $\Omega_{aa}$ couplings. This is, however, constrained by the fact that they would contribute to the EDM of the electron through Barr-Zee diagrams, which can easily be computed by adapting results from the literature. First, the Barr-Zee diagram involving the photon and the Higgs boson of Fig. 2(a) can be computed by adapting the results of Ref. Nakai:2016atk . This diagram and its variations lead to a contribution of

	$\displaystyle\frac{d^{Ah}_{e}}{e}=$	$\displaystyle-\sum_{a}\frac{\alpha\tilde{Q}^{2}_{aa}m_{a}m_{e}}{16\pi^{3}m_{h}^{2}v}\text{Im}(\Omega_{aa})$		(27)
		$\displaystyle\times\int_{0}^{1}dx\frac{1}{x(1-x)}j\left(0,\frac{m_{a}^{2}}{x(1-x)m_{h}^{2}}\right),$

where

$$j(r,s)=\frac{1}{r-s}\left(\frac{r\ln r}{r-1}-\frac{s\ln s}{s-1}\right).$$

(28)

Second, the contribution from the diagrams involving a $Z$ boson and a Higgs boson like that of Fig. 2(b) can also be computed using Ref. Nakai:2016atk and lead to

	$\displaystyle\frac{d^{Zh}_{e}}{e}=\sum_{a,b}\frac{\tilde{Q}_{bb}}{32\pi^{4}m_{h}^{2}}g_{ee}^{V}g_{ee}^{S}\Bigl{(}$	$\displaystyle m_{a}C_{ab}^{1}f_{1}(m_{a},m_{b})$		(29)
	$\displaystyle+$	$\displaystyle m_{b}C_{ab}^{2}f_{2}(m_{a},m_{b})\Bigr{)},$

where

	$\displaystyle C_{ab}^{1}$	$\displaystyle=\text{Re}\left(ig_{ba}^{S}g_{ab}^{A}-g_{ba}^{P}g_{ab}^{V}\right),$		(30)
	$\displaystyle C_{ab}^{2}$	$\displaystyle=-\text{Re}\left(ig_{ba}^{S}g_{ab}^{A}+g_{ba}^{P}g_{ab}^{V}\right),$

with

	$\displaystyle g_{ee}^{S}$	$\displaystyle=-\frac{m_{e}}{v},\quad\quad g_{ee}^{V}=-\frac{\sqrt{g^{2}+{g^{\prime}}^{2}}}{2}\left(-\frac{1}{2}+2s_{W}^{2}\right),$		(31)
	$\displaystyle g_{ab}^{S}$	$\displaystyle=-\frac{(\Omega_{ba}+\Omega_{ab}^{\ast})}{2},\quad g_{ab}^{P}=-i\frac{(\Omega_{ba}-\Omega_{ab}^{\ast})}{2},$
	$\displaystyle g_{ab}^{V}$	$\displaystyle=-\frac{\sqrt{g^{2}+{g^{\prime}}^{2}}}{2}\left(B_{Rba}+B_{Lba}\right),$
	$\displaystyle g_{ab}^{A}$	$\displaystyle=-\frac{\sqrt{g^{2}+{g^{\prime}}^{2}}}{2}\left(B_{Rba}-B_{Lba}\right),$

and

	$\displaystyle f_{1}(m_{a},m_{b})$	$\displaystyle=\int_{0}^{1}dxj\left(\frac{m_{Z}^{2}}{m_{h}^{2}},\frac{\tilde{\Delta}_{ab}}{m_{h}^{2}}\right),$		(32)
	$\displaystyle f_{2}(m_{a},m_{b})$	$\displaystyle=\int_{0}^{1}dxj\left(\frac{m_{Z}^{2}}{m_{h}^{2}},\frac{\tilde{\Delta}_{ab}}{m_{h}^{2}}\right)\frac{1-x}{x},$

with

$$\tilde{\Delta}_{ab}=\frac{xm_{a}^{2}+(1-x)m_{b}^{2}}{x(1-x)}.$$

(33)

Third, the contribution from the diagrams involving two $W$ bosons like in Fig. 2(c) can be computed by adapting the results of Ref. Chang:2005ac . The only subtlety is that the fermions can potentially flow in both directions in the fermion loop. This was not the case in Ref. Chang:2005ac . Thankfully, both diagrams can easily be related and the sum of the two diagrams is

		$\displaystyle\frac{d_{e}^{WW}}{e}=-\frac{\alpha^{2}m_{e}}{8\pi^{2}s_{W}^{4}m_{W}^{2}}\sum_{a,b}\frac{m_{a}m_{b}}{m_{W}^{2}}\text{Im}\left(\hat{A}_{Lba}\hat{A}_{Rba}^{\ast}\right)$		(34)
		$\displaystyle\hskip 42.67912pt\times\left[\tilde{Q}_{bb}\mathcal{G}(r_{a},r_{b},0)+\tilde{Q}_{aa}\mathcal{G}(r_{b},r_{a},0)\right],$

where $r_{a}=m_{a}^{2}/m_{W}^{2}$, $r_{b}=m_{b}^{2}/m_{W}^{2}$, and

		$\displaystyle\mathcal{G}(r_{a},r_{b},r_{c})=$		(35)
		$\displaystyle\quad\int_{0}^{1}\frac{d\gamma}{\gamma}\int_{0}^{1}dyy\Biggl{[}\frac{(R-3K_{ab})R+2(K_{ab}+R)y}{4R(K_{ab}-R)^{2}}$
		$\displaystyle\hskip 76.82234pt+\frac{K_{ab}(K_{ab}-2y)}{2(K_{ab}-R)^{3}}\ln\frac{K_{ab}}{R}\Biggr{]},$

with

$$R=y+(1-y)r_{c},\qquad K_{ab}=\frac{r_{a}}{1-\gamma}+\frac{r_{b}}{\gamma}.$$

(36)

The total EDM is then the sum of Eqs. (27), (29), and (34). In practice, the contribution from Eq. (27) generally overwhelmingly dominates when $\text{BR}(h\to AA^{\prime})$ is close to its maximally allowed value. In this case, the diagram of Fig. 2(b) is suppressed by the fact that $g_{ee}^{V}$ is small due to an accidental partial cancellation and the diagram of Fig. 2(c) is suppressed because the Yukawa couplings $A_{L}$ and $A_{R}$ are simply much larger than $g$. The upper limit on the electron EDM that we use is $|d_{e}|<4.1\times 10^{-30}e\;\text{cm}$ at 90% CL Roussy:2022cmp .⁴⁴4This limit is updated with respect to Ref. PhysRevLett.130.141801 which used Ref. ACME:2018yjb , as the result of Ref. Roussy:2022cmp was not available yet. The impact of this new result on our limits is negligible.

Finally, when $\text{BR}(h\to AA^{\prime})$ is close to the upper limit that we find, the electron EDM generally constrains $A_{L}$ and $A_{R}$ to have a phase difference very close to 0 (or $\pi$). This makes it essentially redundant to consider any other CP-violating observable.

A sizable $\text{BR}(h\to AA^{\prime})$ requires some of the Yukawa couplings $A_{L}$ and $A_{R}$ to be large. These couplings however have the side effect of causing mixing between fields that are part of different representations of the electroweak gauge groups. This means that a large $\text{BR}(h\to AA^{\prime})$ is at risk of generating large contributions to the oblique parameters Peskin:1990zt . The parameters $S$ and $T$ are computed by using the general results for fermions of Refs. Anastasiou:2009rv ; Lavoura:1992np ; Chen:2003fm ; Carena:2007ua ; Chen:2017hak ; Cheung:2020vqm . They are given by:

	$\displaystyle S=\frac{1}{2\pi}\sum_{a,b}\Bigl{\{}$	$\displaystyle\Bigl{(}|\hat{A}_{Lab}|^{2}+|\hat{A}_{Rab}|^{2}\Bigr{)}\psi_{+}(y_{a},y_{b})$		(37)
		$\displaystyle+2\text{Re}\Bigl{(}\hat{A}_{Lab}\hat{A}^{\ast}_{Rab}\Bigr{)}\psi_{-}(y_{a},y_{b})$
		$\displaystyle-\frac{1}{2}\Bigl{[}\Bigl{(}|X_{ab}|^{2}+|X_{Rab}|^{2}\Bigr{)}\chi_{+}(y_{a},y_{b})$
		$\displaystyle+2\text{Re}\Bigl{(}X_{Lab}X^{\ast}_{Rab}\Bigr{)}\chi_{-}(y_{a},y_{b})\Bigr{]}\Bigr{\}},$
	$\displaystyle T=\frac{1}{16\pi s_{W}^{2}c_{W}^{2}}$	$\displaystyle\sum_{a,b}\Bigl{\{}\Bigl{(}|\hat{A}_{Lab}|^{2}+|\hat{A}_{Rab}|^{2}\Bigr{)}\theta_{+}(y_{a},y_{b})$
		$\displaystyle+2\text{Re}\Bigr{(}\hat{A}_{Lab}\hat{A}^{\ast}_{Rab}\Bigr{)}\theta_{-}(y_{a},y_{b})$
		$\displaystyle-\frac{1}{2}\Bigl{[}\Bigl{(}|X_{Lab}|^{2}+|X_{Rab}|^{2}\Bigr{)}\theta_{+}(y_{a},y_{b})$
		$\displaystyle+2\text{Re}\Bigl{(}X_{Lab}X^{\ast}_{Rab}\Bigr{)}\theta_{-}(y_{a},y_{b})\Bigr{]}\Bigr{\}},$

where $y_{a}=m_{a}^{2}/m_{Z}^{2}$, $X_{L/R}=-2B_{L/R}+2\tilde{Q}s_{W}^{2}$ and

	$\displaystyle\psi_{+}(y_{1},y_{2})$	$\displaystyle=\frac{1}{3}-\frac{1}{9}\ln\frac{y_{1}}{y_{2}},$		(38)
	$\displaystyle\psi_{-}(y_{1},y_{2})$	$\displaystyle=-\frac{y_{1}+y_{2}}{6\sqrt{y_{1}y_{2}}},$
	$\displaystyle\chi_{+}(y_{1},y_{2})$	$\displaystyle=\frac{5(y_{1}^{2}+y_{2}^{2})-22y_{1}y_{2}}{9(y_{1}-y_{2})^{2}}$
		$\displaystyle+\frac{3y_{1}y_{2}(y_{1}+y_{2})-y_{1}^{3}-y_{2}^{3}}{3(y_{1}-y_{2})^{3}}\ln\frac{y_{1}}{y_{2}},$
	$\displaystyle\chi_{-}(y_{1},y_{2})$	$\displaystyle=-\sqrt{y_{1}y_{2}}\Bigl{[}\frac{y_{1}+y_{2}}{6y_{1}y_{2}}-\frac{y_{1}+y_{2}}{(y_{1}-y_{2})^{2}}$
		$\displaystyle\hskip 51.21504pt+\frac{2y_{1}y_{2}}{(y_{1}-y_{2})^{3}}\ln\frac{y_{1}}{y_{2}}\Bigr{]},$
	$\displaystyle\theta_{+}(y_{1},y_{2})$	$\displaystyle=y_{1}+y_{2}-\frac{2y_{1}y_{2}}{y_{1}-y_{2}}\ln\frac{y_{1}}{y_{2}},$
	$\displaystyle\theta_{-}(y_{1},y_{2})$	$\displaystyle=2\sqrt{y_{1}y_{2}}\left[\frac{y_{1}+y_{2}}{y_{1}-y_{2}}\ln\frac{y_{1}}{y_{2}}-2\right].$

We use the measurements of the oblique parameters of Ref. Zyla:2020zbs given by

$$S=0.00\pm 0.07,\qquad T=0.05\pm 0.06,$$

(39)

with a correlation of 0.92. We keep points whose $\chi^{2}$ differ by less than 5.99 from the best fit, which corresponds to 95% CL limits.

As a final constraint, the parameters $A_{R}$ and $A_{L}$ are bounded by unitarity. Consider a given scattering between mediators via Higgs exchange and its amplitude $\mathcal{M}$. The latter can be expanded in partial waves as

$$\mathcal{M}=16\pi\sum_{\ell}(2\ell+1)a_{\ell}P_{\ell}(\cos\theta),$$

(40)

where $P_{\ell}(\cos\theta)$ are the Legendre polynomials.