# Neutrino Phenomenology in a Model with Generalized CP symmetry within Type-I
seesaw framework

Tapender<EMAIL_ADDRESS>Department of Physics and Astronomical Science,
Central University of Himachal Pradesh, Dharamshala-176215, INDIA  Sanjeev
Kumar<EMAIL_ADDRESS>Department of Physics and Astrophysics,
University of Delhi, Delhi-110007, INDIA  Surender Verma<EMAIL_ADDRESS>Corresponding Author Department of Physics and Astronomical Science, Central
University of Himachal Pradesh, Dharamshala-176215, INDIA

###### Abstract

We investigate the consequences of generalized CP (GCP) symmetry within the
context of the two Higgs doublet model (2HDM), specifically focusing on the
lepton sector. Utilizing the Type-I seesaw framework, we study an intriguing
connection between the Dirac Yukawa couplings originating from both Higgs
fields, leading to a reduction in the number of independent Yukawa couplings
and simplifying the scalar and Yukawa sectors when compared to the general
2HDM. The CP3 constraint results in two right-handed neutrinos having equal
masses and leads to a diagonal right-handed Majorana neutrino mass matrix.
Notably, CP symmetry experiences a soft break due to the phase associated with
the vacuum expectation value of the second Higgs doublet. The model aligns
well with observed charged lepton masses and neutrino oscillation data,
explaining both masses and mixing angles, and yields distinct predictions for
normal and inverted neutrino mass hierarchies. It features a novel interplay
between atmospheric mixing angle $\theta_{23}$ and neutrino mass hierarchy:
the angle $\theta_{23}$ is below maximal for the normal hierarchy and above
maximal for inverted hierarchy. Another interesting feature of the model is
inherent CP violation for the inverted hierarchy.

##  1 Introduction

The standard model (SM) of particle physics provides a unified and well-tested
theoretical framework for explaining the interactions of known fundamental
particles. It explains how quarks and charged leptons acquire mass. However,
it cannot account for the non-zero mass of neutrinos, which is necessary to
explain observed neutrino oscillations. One way to naturally explain non-zero
neutrino masses is by introducing right-handed neutrinos into the particle
content of the SM and allowing them to have a Majorana mass term. This is
commonly known as the Type-I seesaw mechanism [1, 2, 3, 4, 5]. The smallness
of the neutrino masses can be attained by setting Majorana neutrino mass at
high energy scale.

Extending beyond the Standard Model (SM), a natural step involves adding
another Higgs doublet, known as the Two-Higgs Doublet Model (2HDM). Initially
proposed to address matter-antimatter asymmetry alongside the quark mixing
matrix [6], the 2HDM doesn’t explain neutrino mass. The vacuum expectation
values of these two SU(2) doublets, spontaneously break the CP symmetry
contributing as an extra source for generating matter-antimatter asymmetry[6].
Further, the need for a second Higgs doublet arises naturally in the Minimal
Supersymmetric Standard Model (MSSM) [7] and axion models [8, 9]. Another
reason for considering 2HDM is that it preserves the $\rho$ parameter [6],
connecting the mechanism of electroweak symmetry breaking with the masses of
SM gauge bosons [10].

Despite these characteristic features, 2HDMs have shortcomings, including the
inability to explain neutrino mass, dark matter, and the allowance of tree-
level flavor changing neutral currents (FCNC). The presence of FCNC arises
because both SU(2) scalar doublets can couple to fermions. However, there are
studies suggesting mechanisms to mitigate FCNC interactions. For example:

  1. 1.

FCNC interactions can be fine-tuned by carefully selecting Yukawa couplings
that are suppressed by the heavy mass of the scalar boson responsible for FCNC
[11].

  2. 2.

FCNC can be eliminated by employing a global symmetry, such as Z2, which
restricts a given scalar boson from coupling to fermions of different electric
charges [12, 13].

  3. 3.

Tree-level FCNC can be eliminated by using a global U(1) Peccei-Quinn symmetry
[14].

In addition to addressing FCNC, there have been various attempts to enhance
2HDMs to incorporate neutrino masses [15, 16, 17, 18, 19, 20] and dark matter
[21, 22, 23, 24, 25, 26].

However, 2HDM poses a challenge due to its large number of free parameters,
making it difficult to probe through collider experiments like the LHC. In
general, the scalar potential of the 2HDM consists of fourteen parameters and
can exhibit CP conserving or CP violating behavior [11, 27]. Consequently,
additional constraints are necessary, often derived from symmetry arguments,
to establish relationships among these parameters.

The study of generalized CP (GCP) transformations within the scalar sector of
2HDM is an example of imposing additional symmetries [28, 29]. GCP
transformations can be categorized in various ways. In Ref.[28], they are
classified into three categories: CP1, CP2, and CP3. CP1 and CP2 correspond to
discrete transformations, while CP3 is a continuous transformation that can be
extended to the fermionic sector [29], and they have applied this to the quark
sector. Furthermore, the CP symmetries of the scalar sector in the 2HDM have
been thoroughly investigated using the basis invariant bilinear formalism [30,
31, 32, 33, 34, 35] in many works [31, 32, 33, 34].

In 2HDM extensions addressing neutrino mass, explicit neutrino phenomenology,
including mixing angles and mass-squared differences, is often absent. In this
study, we’ve extended CP3 to the 2HDM’s leptonic Yukawa sector, introducing CP
violation through the second Higgs’s vev phase. Neutrino masses are generated
via the Type-I seesaw relation, involving right-handed neutrinos.

The paper is structured as follows: In Section 2, we present the basic
formalism of 2HDM. Section 3 elaborates on extending CP3 to the neutrino
Yukawa sector within the Type-I seesaw mechanism. We discuss our numerical
analysis in Section 4. Finally, in Section 5, we summarize our conclusions.

##  2 Two Higgs Doublet Model under Generalized CP Symmetry

In 2HDMs, the Standard Model’s field content expands with the addition of an
extra Higgs doublet, denoted as $\Phi_{2}$, which possesses the same charge
assignments as the Higgs field in the SM. However, this minimal expansion in
the scalar sector results in an increased number of free parameters. To
address this parameter growth, it becomes imperative to introduce specific
symmetries. In this context, the Generalized CP (GCP) symmetry is considered.
Under GCP, scalar doublets undergo transformations as elucidated in [28]:

$\Phi_{a}\rightarrow\Phi_{a}^{GCP}=X_{a\alpha}\Phi_{\alpha}^{*}\,\,,$ (1)

where X is an arbitrary unitary CP transformation matrix.

There always exist a choice of basis for which most general GCP transformation
matrix can be brought to the form [36]

$X=\begin{pmatrix}\cos\theta&&\sin\theta\\\
-\sin\theta&&\cos\theta\end{pmatrix}\,\,,$ (2)

where $0\leq\theta\leq\pi/2$. So, we have three distinct cases with respect to
the parameter $\theta$ as mentioned in [28]:

  1. 1.

When $\theta=0$, the symmetry is referred to as CP symmetry of order one
(CP1).

  2. 2.

For $\theta=\pi/2$, the symmetry is known as CP symmetry of order two (CP2).

  3. 3.

In the range $0<\theta<\pi/2$, the symmetry is labeled as CP symmetry of order
three (CP3) and importantly it constitutes a continuous symmetry.

The most general scalar potential with two Higgs doublets can be written as

$\begin{split}V_{H}=&\,\,m^{2}_{11}\Phi_{1}^{\dagger}\Phi_{1}+m^{2}_{22}\Phi_{2}^{\dagger}\Phi_{2}-[m^{2}_{12}\Phi_{1}^{\dagger}\Phi_{2}+H.c.]+\frac{1}{2}\lambda_{1}(\Phi_{1}^{\dagger}\Phi_{1})^{2}+\frac{1}{2}\lambda_{2}(\Phi_{2}^{\dagger}\Phi_{2})^{2}\\\
&+\lambda_{3}(\Phi_{1}^{\dagger}\Phi_{1})(\Phi_{2}^{\dagger}\Phi_{2})+\lambda_{4}(\Phi_{1}^{\dagger}\Phi_{2})(\Phi_{2}^{\dagger}\Phi_{1})\\\
&+\big{[}\frac{1}{2}\lambda_{5}(\Phi_{1}^{\dagger}\Phi_{2})^{2}+\lambda_{6}(\Phi_{1}^{\dagger}\Phi_{1})(\Phi_{1}^{\dagger}\Phi_{2})+\lambda_{7}(\Phi_{2}^{\dagger}\Phi_{2})(\Phi_{1}^{\dagger}\Phi_{2})+H.c.\big{]}\,\,,\end{split}$
(3)

having total 14 parameters. Here, $m_{11}^{2}$, $m_{22}^{2}$ and $\lambda_{1}$
through $\lambda_{4}$ are real parameters while $m_{12}^{2}$, $\lambda_{5}$,
$\lambda_{6}$ and $\lambda_{7}$ are generally complex.

Based on the findings in Ref. [29], it’s established that, in addition to the
standard CP symmetry CP1, CP3 is the only symmetry that can be extended to the
Yukawa sector for leptons. To maintain CP3 invariance within the scalar
potential, certain conditions must be satisfied. Specifically, we must have
$m_{11}^{2}=m_{22}^{2}$, $m_{12}^{2}=0$, $\lambda_{2}=\lambda_{1}$,
$\lambda_{6}=0$, $\lambda_{7}=0$ and
$\lambda_{5}=\lambda_{1}-\lambda_{3}-\lambda_{4}$ which must be a real
parameter.

To avoid the presence of Goldstone bosons following spontaneous symmetry
breaking, it’s necessary to introduce soft CP3 symmetry breaking. Therefore,
we will consider $m_{11}^{2}\neq m_{22}^{2}$ and $\Re[m_{12}^{2}]\neq 0$. This
softly broken CP3 symmetric potential also leads to a CP-violating vacuum
expectation value (vev) for the second Higgs doublet, which assumes a crucial
role in the exploration of CP violation in the lepton sector, as we will delve
into in the upcoming sections.

##  3 CP3 in Yukawa Sector with Type-I seesaw

In our study, we’ve expanded upon the Standard Model (SM) by incorporating one
Higgs doublet and three right-handed neutrinos denoted as $N_{R}$. Within the
framework of the Type-I seesaw mechanism, the relevant Yukawa Lagrangian
responsible for generating the masses of both charged leptons and
neutrinos111For quark masses, see Ref. [29] is expressed as follows

$-\mathcal{L}_{Y}=\overline{L}_{L}\Gamma_{a}\Phi_{a}l_{R}+\overline{L}_{L}Y_{a}\tilde{\Phi}_{a}N_{R}+\frac{1}{2}\overline{N^{c}_{R}}MN_{R}+H.c.\,\,,$
(4)

where $L_{L}$, $l_{R}$ are Standard Model $SU(2)$ left-handed doublets and
right-handed singlets, $\Phi_{a}$ (a = 1, 2) are Higgs doublets and NR are
right-handed neutrino singlets. $\Gamma_{a}$ and $Y_{a}$ are the Yukawa
coupling matrices for charged leptons and neutrinos respectively and M is
lepton number violating Majorana mass term for right-handed neutrinos. Now we
will extend the GCP symmetry to the leptonic Yukawa sector.

The fields involved in Eqn.(4) transforms under GCP symmetry as

$\begin{rcases}\begin{aligned}
\Phi_{a}&\rightarrow\Phi_{a}^{GCP}=X_{ab}\Phi_{b}^{*}\,\,,\\\
\tilde{\Phi}_{a}&\rightarrow\tilde{\Phi}_{a}^{GCP}=X_{ab}^{*}(\tilde{\Phi}_{b}^{\dagger})^{T}\,\,,\\\
L_{L}&\rightarrow
L_{L}^{GCP}=iX_{\zeta}\gamma^{0}C\overline{L}_{L}^{T}\,\,,\\\
l_{R}&\rightarrow
l_{R}^{GCP}=iX_{\beta}\gamma^{0}C\overline{l}_{R}^{T}\,\,,\\\
N_{R}&\rightarrow
N_{R}^{GCP}=iX_{\gamma}\gamma^{0}C\overline{N}_{R}^{T}\,\,,\end{aligned}\end{rcases}$
(5)

where $\gamma^{0}$ is Dirac gamma matrix and C is charge conjugation matrix,
$X$, $X_{\zeta}$, $X_{\beta}$ and $X_{\gamma}$ are CP transformation matrices.

For Lagrangian to remain invariant under these CP transformations we find
Yukawa coupling matrices to transform as

$\Gamma_{b}^{*}=X_{\zeta}^{\dagger}\Gamma_{a}X_{\beta}X_{ab}\,\,,$ (6)
$Y_{b}^{*}=X_{\zeta}^{\dagger}Y_{a}X_{\gamma}X_{ab}^{*}\,\,,$ (7)

and Majorana mass matrix to transform as

$M^{*}=X_{\gamma}^{T}MX_{\gamma}\,\,,$ (8)

where

$M=\begin{pmatrix}M_{11}&&M_{12}&&M_{13}\\\ M_{12}&&M_{22}&&M_{23}\\\
M_{13}&&M_{23}&&M_{33}\end{pmatrix}\,\,.$ (9)

The CP transformation matrices involved are given by

$X_{\zeta}=\begin{pmatrix}\cos\zeta&&\sin\zeta&&0\\\
-\sin\zeta&&\cos\zeta&&0\\\ 0&&0&&1\end{pmatrix}\,\,,$ (10)
$X_{\beta}=\begin{pmatrix}\cos\beta&&\sin\beta&&0\\\
-\sin\beta&&\cos\beta&&0\\\ 0&&0&&1\end{pmatrix}\,\,,$ (11)
$X_{\gamma}=\begin{pmatrix}\cos\gamma&&\sin\gamma&&0\\\
-\sin\gamma&&\cos\gamma&&0\\\ 0&&0&&1\end{pmatrix}\,\,.$ (12)

It was found in Ref. [29] that CP3 symmetry with $\theta=\pi/3$
($\zeta=\beta=\gamma=\pi/3$) can be extended to Yukawa sector producing
correct quark masses. Under these conditions, the forms of Yukawa coupling
matrices given in Eqns.(6) and (7) become

$\Gamma_{1}=\begin{pmatrix}ia_{11}&&ia_{12}&&a_{13}\\\
ia_{12}&&-ia_{11}&&a_{23}\\\ a_{31}&&a_{32}&&0\end{pmatrix}\,\,,\,\,\
\Gamma_{2}=\begin{pmatrix}ia_{12}&&-ia_{11}&&-a_{23}\\\
-ia_{11}&&-ia_{12}&&a_{13}\\\ -a_{32}&&a_{31}&&0\end{pmatrix}\,\,,$ (13)
$Y_{1}=\begin{pmatrix}ib_{11}&&ib_{12}&&b_{13}\\\ ib_{12}&&-ib_{11}&&b_{23}\\\
b_{31}&&b_{32}&&0\end{pmatrix}\,\,,\,\,\
Y_{2}=\begin{pmatrix}ib_{12}&&-ib_{11}&&-b_{23}\\\
-ib_{11}&&-ib_{12}&&b_{13}\\\ -b_{32}&&b_{31}&&0\end{pmatrix}\,\,,$ (14)

where all $a$’s and $b$’s are real parameters. The choice of $\theta=\pi/3$
alongside $\zeta=\beta=\gamma=\pi/3$ in the leptonic sector stems from the
similarity in GCP transformation properties between quarks and leptonic
fields, as outlined in Eqn.(5).

Now, we need to solve constrains given by Eqn.(8) which can be rewritten as:

$M^{*}-X_{\gamma}^{T}MX_{\gamma}=0\,\,.$ (15)

Using Eqns.(9) and (12) in Eqn.(15), the set of constraints are

$\displaystyle M_{11}^{*}-M_{11}\cos^{2}\gamma-M_{22}\sin^{2}\gamma+M_{12}\sin
2\gamma$ $\displaystyle=$ $\displaystyle 0,$ (16) $\displaystyle
M_{12}^{*}-M_{12}\cos 2\gamma+(-M_{11}+M_{22})\sin\gamma\cos\gamma$
$\displaystyle=$ $\displaystyle 0,$ (17) $\displaystyle
M_{22}^{*}-M_{22}\cos^{2}\gamma-2M_{12}\cos\gamma\sin\gamma-
M_{11}\sin^{2}\gamma$ $\displaystyle=$ $\displaystyle 0,$ (18) $\displaystyle
M_{13}^{*}-M_{13}\cos\gamma+M_{23}\sin\gamma$ $\displaystyle=$ $\displaystyle
0,$ (19) $\displaystyle M_{23}^{*}-M_{23}\cos\gamma-M_{13}\sin\gamma$
$\displaystyle=$ $\displaystyle 0,$ (20) $\displaystyle M_{33}^{*}-M_{33}$
$\displaystyle=$ $\displaystyle 0.$ (21)

In Eqns.(16), (17) and (18), the real part can be separated out as

$\begin{pmatrix}1-\cos^{2}\gamma&\sin 2\gamma&-\sin^{2}\gamma\\\
-\cos\gamma\sin\gamma&1-\cos 2\gamma&\cos\gamma\sin\gamma\\\
-\sin^{2}\gamma&-2\cos\gamma\sin\gamma&1-\cos^{2}\gamma\\\
\end{pmatrix}\begin{pmatrix}\Re[M_{11}]\\\ \Re[M_{12}]\\\ \Re[M_{22}]\\\
\end{pmatrix}=0,$ (22)

and the imaginary part can be separated out as

$\begin{pmatrix}-1-\cos^{2}\gamma&\sin 2\gamma&-\sin^{2}\gamma\\\
-\cos\gamma\sin\gamma&-1-\cos 2\gamma&\cos\gamma\sin\gamma\\\
-\sin^{2}\gamma&-2\cos\gamma\sin\gamma&-1-\cos^{2}\gamma\\\
\end{pmatrix}\begin{pmatrix}\Im[M_{11}]\\\ \Im[M_{12}]\\\ \Im[M_{22}]\\\
\end{pmatrix}=0.$ (23)

. Further, from Eqns.(19) and (20) we have, for real part

$\begin{pmatrix}1-\cos\gamma&\sin\gamma\\\ -\sin\gamma&1-\cos\gamma\\\
\end{pmatrix}\begin{pmatrix}\Re[M_{13}]\\\ \Re[M_{23}]\end{pmatrix}=0,$ (24)

and, for imaginary part,

$\begin{pmatrix}-1-\cos\gamma&\sin\gamma\\\ -\sin\gamma&-1-\cos\gamma\\\
\end{pmatrix}\begin{pmatrix}\Im[M_{13}]\\\ \Im[M_{23}]\end{pmatrix}=0.$ (25)

For $\gamma=\pi/3$, we have

$\begin{pmatrix}\frac{3}{4}&\frac{\sqrt{3}}{2}&\frac{-3}{4}\\\
-\frac{\sqrt{3}}{4}&\frac{3}{2}&\frac{\sqrt{3}}{4}\\\
\frac{-3}{4}&\frac{-\sqrt{3}}{2}&\frac{3}{4}\end{pmatrix}\begin{pmatrix}\Re[M_{11}]\\\
\Re[M_{12}]\\\ \Re[M_{22}]\\\ \end{pmatrix}=0,$ (26)
$\begin{pmatrix}\frac{-5}{4}&\frac{\sqrt{3}}{2}&\frac{-3}{4}\\\
-\frac{\sqrt{3}}{4}&-\frac{1}{2}&\frac{\sqrt{3}}{4}\\\
\frac{-3}{4}&-\frac{\sqrt{3}}{2}&-\frac{5}{4}\end{pmatrix}\begin{pmatrix}\Im[M_{11}]\\\
\Im[M_{12}]\\\ \Im[M_{22}]\\\ \end{pmatrix}=0,$ (27)

and

$\begin{pmatrix}\frac{1}{2}&\frac{\sqrt{3}}{2}\\\
-\frac{\sqrt{3}}{2}&\frac{1}{2}\\\ \end{pmatrix}\begin{pmatrix}\Re[M_{13}]\\\
\Re[M_{23}]\end{pmatrix}=0,$ (28)
$\begin{pmatrix}-\frac{3}{2}&\frac{\sqrt{3}}{2}\\\
-\frac{\sqrt{3}}{2}&-\frac{3}{2}\\\ \end{pmatrix}\begin{pmatrix}\Im[M_{13}]\\\
\Im[M_{23}]\end{pmatrix}=0.$ (29)

In Eqns.(27), (28), and (29), the determinant of the square matrix is non-
zero, implying a unique solution where $\Im[M_{11}]$, $\Im[M_{12}]$,
$\Im[M_{22}]$, $\Re[M_{13}]$, $\Re[M_{23}]$, $\Im[M_{13}]$, and $\Im[M_{23}]$
all equal zero. On the other hand, in Eqn.(26), the determinant of the square
matrix is zero, indicating arbitrary solutions, with $\Re[M_{12}]$ equal to
zero and $\Re[M_{11}]\equiv M_{1}$ equal to $\Re[M_{22}]$. Furthermore,
Eqn.(21) leads to $\Im[M_{33}]$ being zero, leaving only $\Re[M_{33}]\equiv
M_{3}$ as the relevant parameter.

It’s worth noting that the CP3 constraint results in two right-handed
neutrinos having equal masses, leading to a diagonal matrix $M$, described as
follows

$M=\begin{pmatrix}M_{1}&&0&&0\\\ 0&&M_{1}&&0\\\ 0&&0&&M_{3}\end{pmatrix}\,\,.$
(30)

After spontaneous symmetry breaking (SSB), both Higgs doublets get vevs, given
by

$\left<\Phi_{1}\right>=\begin{pmatrix}0\\\
\frac{v_{1}}{\sqrt{2}}\end{pmatrix}\,\,\,,\,\,\left<\Phi_{2}\right>=\begin{pmatrix}0\\\
e^{i\alpha}\frac{v_{2}}{\sqrt{2}}\end{pmatrix}\,\,,$ (31)

with condition that $v=\sqrt{v_{1}^{2}+v_{2}^{2}}\approx 245$ GeV, where $v$
is the standard model vev. Consequently, the charged leptons mass matrix
becomes

$\displaystyle M_{l}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(v_{1}\Gamma_{1}+e^{i\alpha}v_{2}\Gamma_{2}),$
(32) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(\cos\phi\Gamma_{1}+e^{i\alpha}\sin\phi\Gamma_{2})v,$
(33)

and for neutrinos we have, Dirac mass matrix given by

$\displaystyle M_{D}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(v_{1}Y_{1}+e^{-i\alpha}v_{2}Y_{2}),$ (34)
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2}}(\cos\phi
Y_{1}+e^{-i\alpha}\sin\phi Y_{2})v,$ (35)

where $\phi$ is defined as $\tan\phi=v_{2}/v_{1}$.

We work in the basis in which charged lepton mass matrix is diagonal. The
charged lepton mass matrix can be diagonalized as

$M_{l}=U_{l}m_{diag}U_{R}^{\dagger},$ (36)

where $U_{l}$ and $U_{R}$ are $3\times 3$ unitary matrices and
$m_{diag}=diag(m_{e},m_{\mu},m_{\tau})$ is diagonal matrix with positive real
entries giving mass eigenvalues of electron, muon, and tau, respectively. So,
we have

$\displaystyle U_{l}^{\dagger}M_{l}M_{l}^{\dagger}U_{l}=m_{diag}^{2},$ (37)

where $U_{l}$ rotates $M_{D}$ into the basis in which charged leptons are
diagonal

$M_{D}^{new}=U_{l}^{\dagger}M_{D}.$ (38)

Using Type-I seesaw, the effective light neutrinos mass matrix is given by

$\displaystyle M_{\nu}$ $\displaystyle=$ $\displaystyle-
M_{D}^{new}M^{-1}(M_{D}^{new})^{T}\,,$ (39) $\displaystyle=$
$\displaystyle-(U_{l}^{\dagger}M_{D})M^{-1}(U_{l}^{\dagger}M_{D})^{T},$ (40)

which is a complex symmetric matrix. This matrix is related to Yukawa coupling
matrices $Y_{1}$ and $Y_{2}$ through Eqn.(35). The effective light neutrino
mass matrix can be diagonalized by $3\times 3$ unitary matrix $U$ as

$U^{\dagger}M_{\nu}U^{*}=m,$ (41)

where $m_{ik}=m_{i}\delta_{ik}$, $m_{i}>0$ $(i,k=1,2,3)$.

We will now move forward with the numerical determination of charged lepton
masses and the parameters governing neutrino oscillations. This process
entails the variation of free parameters to ascertain the permissible
parameter space within the model.

##  4 Numerical Analysis and Discussion

In our numerical analysis, we generated random numbers uniformly for the vev-
phase $\alpha$ in the range of 0 to 2$\pi$. We also generated random numbers
uniformly for the masses of the right-handed neutrinos $M_{1}$ and $M_{3}$,
which ranged from 1011 to 1013 GeV and from 1.1$\times$1013 to 1015 GeV,
respectively. We considered two cases for the vev $v_{1}$, as discussed in the
following subsections.

###  4.1 When $v_{1}<<v_{2}$

In this scenario, we examined the influence of a very small vev on our
parameter space. To emphasize the dominance of vev $v_{2}$, we randomly varied
$v_{1}$ in the range of $(0-5)\times 10^{-6}$ GeV. The value of $v_{2}$ is
subsequently determined using the equation $v_{2}=\sqrt{v^{2}-v_{1}^{2}}$ GeV.
We then determined the masses of charged leptons and the parameters governing
neutrino oscillations, which we discuss in the following subsections.

####  4.1.1 Charged Lepton Masses

To compute the masses of charged leptons, we varied the Yukawa coupling
parameters within the range specified in Table 1. We then proceeded to
numerically diagonalize the mass matrix $M_{l}M_{l}^{\dagger}$, as described
in Eqn.(37), to obtain the squared masses of charged leptons
($m_{e}^{2},m_{\mu}^{2},m_{\tau}^{2}$). Through this analysis, we identified
parameter values that consistently yielded the correct charged lepton masses
for both normal ($m_{1}<m_{2}<m_{3}$) and inverted ($m_{3}<m_{1}<m_{2}$)
hierarchies of neutrinos. The benchmark points are listed in the second column
of Table 2. With these parameter values, we calculated the charged lepton
masses, as presented in Table 3, which closely align with experimentally
observed values.

Yukawa Coupling |  When $v_{1}<<v_{2}$ | When $v_{1}$ ranges from ($10-17$) Gev
---|---|---
Parameter |  Range |  Range
$a_{11}$ | $3\times 10^{-5}-6\times 10^{-5}$ |  $2\times 10^{-5}-3\times 10^{-5}$
$a_{12}$ | $8\times 10^{-5}-1\times 10^{-4}$ |  $6\times 10^{-5}-8\times 10^{-5}$
$a_{13}$ | $6\times 10^{-3}-8\times 10^{-3}$ |  $6\times 10^{-3}-7\times 10^{-3}$
$a_{23}$ | $6\times 10^{-3}-8\times 10^{-3}$ |  $8\times 10^{-3}-9\times 10^{-3}$
$a_{31}$ | $4\times 10^{-4}-6\times 10^{-4}$ |  $4\times 10^{-4}-6\times 10^{-4}$
$a_{32}$ | $1\times 10^{-4}-2\times 10^{-4}$ |  $1\times 10^{-4}-2\times 10^{-4}$
Table 1: The ranges of Yukawa couplings used in numerical analysis for charged leptons. |
---|---
Figure 1: Predictions for effective Majorana mass $|m_{ee}|$ for normal (left) and inverted (right) hierarchy when $v_{1}<<v_{2}$. |
---|---
Figure 2: $\delta-J_{CP}$ correlation with vev-phase $\alpha=0$ and $\alpha\neq 0$ for normal (left) and inverted (right) hierarchies. |
---|---
Figure 3: Correlation between the mixing angles $\theta_{13}$ and $\theta_{23}$ for normal (left) and inverted (right) hierarchies. |
---|---
|
Figure 4: Correlation plots between the effective neutrino mass $|m_{ee}|$ and lightest neutrino mass, as well as between $J_{CP}$ and CP-violating phase $\delta$, for both normal (first row) and inverted (second row) hierarchies when the vev $v_{1}$ ranges from 10 to 17 GeV. |
---|---
Figure 5: Correlation between the mixing angles $\theta_{13}$ and $\theta_{23}$ for normal (left) and inverted (right) hierarchies. |
---|---
|
Figure 6: Correlation of $\delta$ with $\theta_{23}$ and vev-phase $\alpha$ for NH (first row) and IH (second row). Parameters |  When $v_{1}<<v_{2}$ |  When $v_{1}$ ranges from ($10-17$) GeV
---|---|---
(Units) |  | NH | IH
$a_{11}$ | $5.740882\times 10^{-5}$ | $2.692411\times 10^{-5}$ | $2.692412\times 10^{-5}$
$a_{12}$ | $9.066478\times 10^{-5}$ | $7.993756\times 10^{-5}$ | $7.993757\times 10^{-5}$
$a_{13}$ | $7.124930\times 10^{-3}$ | $6.382725\times 10^{-3}$ | $6.382726\times 10^{-3}$
$a_{23}$ | $7.377109\times 10^{-3}$ | $8.028028\times 10^{-3}$ | $8.028028\times 10^{-3}$
$a_{31}$ | $5.843790\times 10^{-4}$ | $5.945601\times 10^{-4}$ | $5.945602\times 10^{-4}$
$a_{32}$ | $1.377038\times 10^{-4}$ | $1.065540\times 10^{-4}$ | $1.065540\times 10^{-4}$
$v_{1}$ (GeV) | $3.287415\times 10^{-6}$ | $11.00$ | $15.30$
$v_{2}$ (GeV) | $245$ | $244.75$ | $244.52$
$\alpha$ (o) | $197.54$ | $177.68$ | $170.20$
Table 2: Benchmark points for both the cases (i) $v_{1}<<v_{2}$ (second column) (ii) $v_{1}$ ranges from ($10-17$) GeV (third column), yielding correct values of charged lepton masses. Masses |  When $v_{1}<<v_{2}$ |  When $v_{1}$ ranges from ($10-17$) Gev
---|---|---
(MeV) |  | NH | IH
$m_{e}$ | $0.510$ | $0.511$ | $0.511$
$m_{\mu}$ | $105.66$ | $105.66$ | $105.66$
$m_{\tau}$ | $1776.87$ | $1776.85$ | $1776.85$
Table 3:  The values of the charged lepton masses obtained using benchmark points given in Table 2. Yukawa |  When $v_{1}<<v_{2}$ |  When $v_{1}$ ranges from ($10-17$) Gev
---|---|---
Coupling | Range | Range
Parameter |  NH |  IH |  NH |  IH
$b_{11}$ | $6\times 10^{-2}-8\times 10^{-2}$ | $7\times 10^{-2}-2\times 10^{-1}$ | $6\times 10^{-2}-8\times 10^{-2}$ | $8\times 10^{-2}-9\times 10^{-2}$
$b_{12}$ | $1\times 10^{-2}-5\times 10^{-2}$ | $2\times 10^{-3}-4\times 10^{-3}$ | $2\times 10^{-2}-3\times 10^{-2}$ | $2\times 10^{-3}-3\times 10^{-3}$
$b_{13}$ | $7\times 10^{-1}-9\times 10^{-1}$ | $3\times 10^{-1}-5\times 10^{-1}$ | $6\times 10^{-1}-8\times 10^{-1}$ | $1\times 10^{-1}-2\times 10^{-1}$
$b_{23}$ | $1\times 10^{-4}-5\times 10^{-4}$ | $4\times 10^{-1}-6\times 10^{-1}$ | $1\times 10^{-3}-2\times 10^{-3}$ | $4\times 10^{-1}-6\times 10^{-1}$
$b_{31}$ | $1\times 10^{-4}-5\times 10^{-4}$ | $3\times 10^{-2}-5\times 10^{-2}$ | $1\times 10^{-2}-2\times 10^{-2}$ | $3\times 10^{-2}-5\times 10^{-2}$
$b_{32}$ | $4\times 10^{-2}-6\times 10^{-2}$ | $3\times 10^{-2}-5\times 10^{-2}$ | $4\times 10^{-2}-6\times 10^{-2}$ | $6\times 10^{-2}-7\times 10^{-2}$
Table 4:  The ranges of Yukawa couplings used in numerical analysis for normal and inverted hierarchies of neutrinos. Parameter |  When $v_{1}<<v_{2}$ |  When $v_{1}$ ranges from ($10-17$) Gev
---|---|---
(Units) |  NH |  IH |  NH |  IH
$b_{11}$ | $7.244577\times 10^{-2}$ | $1.066963\times 10^{-1}$ | $6.009812\times 10^{-2}$ | $8.006565\times 10^{-2}$
$b_{12}$ | $1.498695\times 10^{-2}$ | $3.644816\times 10^{-3}$ | $2.334813\times 10^{-2}$ | $2.601354\times 10^{-3}$
$b_{13}$ | $8.560231\times 10^{-1}$ | $4.519143\times 10^{-1}$ | $6.935666\times 10^{-1}$ | $1.522807\times 10^{-1}$
$b_{23}$ | $1.033423\times 10^{-4}$ | $5.112375\times 10^{-1}$ | $1.791234\times 10^{-3}$ | $4.531115\times 10^{-1}$
$b_{31}$ | $1.356340\times 10^{-4}$ | $4.546001\times 10^{-2}$ | $1.655557\times 10^{-2}$ | $3.909223\times 10^{-2}$
$b_{32}$ | $5.796197\times 10^{-2}$ | $4.435139\times 10^{-2}$ | $4.880238\times 10^{-2}$ | $6.157254\times 10^{-2}$
$M_{1}$ (GeV) | $2.447917\times 10^{12}$ | $6.721219\times 10^{12}$ | $3.066113\times 10^{12}$ | $4.417490\times 10^{12}$
$M_{3}$ (GeV) | $2.222544\times 10^{14}$ | $5.140671\times 10^{14}$ | $2.748449\times 10^{14}$ | $6.211162\times 10^{14}$
Table 5: The benchmark points yielding correct neutrino phenomenology (i.e. neutrino mixing angles and mass-squared differences are within $3\sigma$ experimental range[37]) for normal and inverted hierarchical neutrino masses. Parameters |  When $v_{1}<<v_{2}$ |  When $v_{1}$ ranges from ($10-17$) Gev
---|---|---
(Units) |  NH |  IH |  NH |  IH
$\theta_{13}$ (o) | $8.44$ | $8.85$ | $8.86$ | $8.22$
$\theta_{12}$ (o) | $31.70$ | $34.30$ | $32.10$ | $31.80$
$\theta_{23}$ (o) | $43.76$ | $47.09$ | $42.96$ | $48.14$
$\Delta m_{21}^{2}$ (eV2) | $7.61\times 10^{-5}$ | $7.33\times 10^{-5}$ | $7.64\times 10^{-5}$ | $7.04\times 10^{-5}$
$\Delta m_{31}^{2}$ (eV2) | $2.58\times 10^{-3}$ | $2.48\times 10^{-3}$ | $2.49\times 10^{-3}$ | $2.43\times 10^{-3}$
Table 6: The values of the neutrino oscillation parameters obtained using
benchmark points given in Table 5 for normal and inverted hierarchical
neutrino masses.

####  4.1.2 Neutrino Oscillation Parameters

To find the neutrino mixing angles and mass-squared differences, we first
rotated the neutrino mass matrix using the previously derived $U_{l}$. The
matrix $M_{\nu}$, described in Eqn.(40), generally contains Yukawa couplings
$Y_{1}$ and $Y_{2}$, as given in Eqn.(14). We then introduced random
variations to these Yukawa coupling parameters within the specified ranges
presented in the second column of Table 4, considering both normal (NH) and
inverted (IH) hierarchies. This process allowed us to identify a benchmark
point that consistently yielded values for the neutrino oscillation
parameters[37]. These values are listed in the second column of Table 5 for
reference. The corresponding values of the mixing angles and mass-squared
differences can be found in the second column of Table 6.

For these parameter values, the effective Majorana neutrino mass takes the
values $|m_{ee}|=\left|\sum_{i}U_{ei}^{2}m_{i}\right|=0.02693$ eV and
$|m_{ee}|=0.04952$ eV for NH and IH of neutrinos, respectively. The masses
$m_{1}$, $m_{2}$, and $m_{3}$ exhibited significant degeneracy in the case of
NH, while for IH, there was an order of magnitude difference between the
lightest mass and the other heavier masses. A linear correlation between
$|m_{ee}|$ and the lightest neutrino mass ($m_{1}$ for NH and $m_{3}$ for IH)
is evident in Fig. 1.

CP violation in the leptonic sector remains unobserved to date. Consequently,
it will be instructive to scrutinize the model’s predictions for CP violation.
In our model, CP violation, effectively, seeds from the complex vev
$\left<\Phi_{2}\right>$ via its phase $\alpha$. CP violation can be understood
in a rephasing invariant way by defining Jarlskog parameter
$J_{CP}=\Im\left(U_{\mu 3}U_{e3}^{*}U_{e2}U_{\mu 2}^{*}\right)$ [38, 39, 40].
For the present case, i.e. $v_{1}<<v_{2}$, we find that $J_{CP}$ and $\delta$
(the phase of the $U_{e3}$ element in $U$) are exceedingly small, regardless
of the neutrino mass hierarchy, and whether the vev-phase $\alpha$ is equal to
zero or non-zero (see Fig. (2)). Figure (2) clearly illustrates that when
$v_{1}<<v_{2}$, the model predicts a effective CP conserving scenario,
regardless of the value of the vev-phase $\alpha$.

Furthermore, in this parameter space region, the mixing angle $\theta_{23}$
falls within the lower octant ($\theta_{23}<45^{o}$) for the normal hierarchy
(NH) and the upper octant ($\theta_{23}>45^{o}$) for the inverted hierarchy
(IH), as illustrated in Fig. 3. Specifically, we note that the mixing angle
$\theta_{23}$ is approximately $44^{o}$ for the normal hierarchy and $47^{o}$
for the inverted hierarchy. In the subsequent section, we will explore an
alternative scenario by increasing the vev $v_{1}$ to observe its impact on CP
violation in the leptonic sector.

###  4.2 When vev $v_{1}$ is in the GeV range

In this specific scenario, we randomly varied the vev $v_{1}$ within the range
of 10 to 17 GeV to analyze the influence of both vevs on the parameter space.
Once again, in this case, $v_{2}$ is determined by the relation
$v_{2}=\sqrt{v^{2}-v_{1}^{2}}$ GeV. The predictions obtained for masses and
other parameters will be discussed as follows.

####  4.2.1 Charged lepton masses

In this case we have randomly varied the Yukawa couplings in range as shown in
third column of Table 1 and after diagonalizing the mass matrix
$M_{l}M_{l}^{\dagger}$ numerically using Eqn.(37) we have obtained squared-
masses of the charged leptons. The benchmark point and prediction for the
masses of charged leptons are listed in third column of Tables 2 and 3,
respectively.

####  4.2.2 Neutrino Oscillation Parameters

The ranges of Yukawa couplings considered in the numerical analysis are
presented in the third column of Table 4. In both normal and inverted
hierarchies of neutrinos, we have identified the benchmark point and obtained
the values of the mixing angles and mass squared differences, as shown in the
third column of Tables 5 and 6, respectively.

For these values, we have determined that $|m_{ee}|$ is approximately 0.01191
eV for NH and 0.04726 eV for IH. Even with the increased value of $v_{1}$, the
masses $m_{1}$, $m_{2}$, and $m_{3}$ follow the same trend as in the previous
case when $v_{1}<<v_{2}$.

The correlation plots, in Fig. 4, are shown for NH (first row) and IH (second
row) cases. It can be seen from Fig. 4 that $|m_{ee}|$ is strongly correlated
to $m_{1}$ (Fig. 4(a)) and Jarlskog invariant $J_{CP}\in(-0.006\rightarrow
0.006)$ (Fig. 4(b)) for NH. For IH $|m_{ee}|\in(4.67\rightarrow 4.83)\times
10^{-2}$ eV (Fig. 4(c)) and
$J_{CP}\in(-0.021\rightarrow-0.018)\oplus(0.018\rightarrow 0.021)$ (Fig.
4(d)). It is interesting to note that in case of IH, the CP violating phase
$\delta=0$ is disallowed thus making this scenario necessarily CP violating.
Another characteristic feature of the model is predicted correlation between
yet unknown $\theta_{23}$-octant and neutrino mass hierarchy. In Fig. 5 we
have shown allowed parameter space of the model in ($\theta_{13}-\theta_{23}$)
plane. It is evident from these plots that, for NH (IH), mixing angles
($\theta_{13}-\theta_{23}$) exhibit negative (positive) correlation. Thus,
$\theta_{23}$ resides in the lower octant for NH and in higher octant for IH
see Fig. 5. Further, the model exhibit a sharp prediction for $\theta_{23}$
approximately equal to $43^{o}$ ($48^{o}$) for NH (IH). Similarly,
$\theta_{13}$ is found to be around $8.8^{o}$ ($8.24^{o}$) for NH (IH). Some
other interesting correlation plots are shown in Fig. 6.

Multi-Higgs models often lead to significant flavor-changing neutral currents
(FCNC)[41]. One approach to reduce these couplings is by aligning all right-
handed fermions to interact with a single Higgs. This alignment can be
achieved through an additional global $Z_{2}$ symmetry[12, 13]. Alternatively,
it can result from adjusting the ratio of the vevs $v_{2}/v_{1}$. In our
study, we are investigating the consequences of Generalized CP (GCP) symmetry
without imposing additional symmetry constraints on the most general Yukawa
couplings and their corresponding neutrino phenomenology. Within the model it
is possible that the Higgs responsible for FCNC is very heavy and thus,
produce vanishing FCNCs. The examination of FCNC effects in the current model
is beyond the scope of this study and will be explored in future research.

##  5 Conclusions

In extended theoretical frameworks beyond the Standard Model, the count of
free parameters increases compared to those at low energy. By introducing
additional symmetry into the Lagrangian, we can substantially reduce the
number of free parameters. In this study, we investigate the impact of GCP
symmetries in the leptonic sector within the context of the 2HDM model. To
incorporate non-zero neutrino masses, we extend the lepton sector by
introducing right-handed neutrinos through the Type-I seesaw mechanism.

The scalar potential of the 2HDM model typically involves fourteen free
parameters. However, due to the GCP symmetry we impose, this number reduces to
four in the unbroken CP3 case and six in the softly broken CP3 symmetry case.
Consequent to GCP, the charged lepton (Eqn. (13)) and neutrino Yukawa coupling
(Eqn. (14)) matrices contain 12 independent parameters, six each in charged
lepton and neutrino sectors. Also, Majorana mass term have two real parameters
(Eqn. (30)).

This model exhibits a rich phenomenology and reveals strong correlations among
neutrino oscillation parameters. The complex vev phase $\alpha$ is the sole
source of CP violation in the model. We consider two distinct phenomenological
scenarios:

  1. 1.

In the first scenario, where vev $v_{1}$ is much smaller than $v_{2}$, CP is
conserved regardless of the value of the vev-phase $\alpha$ (see Fig. 2). This
scenario provides a unique phenomenology for normal and inverted hierarchies
of neutrinos. Notably, the model precisely predicts the neutrino mixing
angles, particularly the atmospheric mixing angle $\theta_{23}$.

  2. 2.

In the second scenario, where $v_{1}$ is in the GeV range, the atmospheric
mixing angle $\theta_{23}$ is below (for NH) or above (for IH) maximality,
approximately $\approx 43^{o}$ and $\approx 48^{o}$ respectively (see Fig. 5).
The Dirac CP phase $\delta$ is tightly constrained to be within the range of
$-10^{o}$ to $10^{o}$ for the NH case. However, if the neutrino masses follow
an inverted mass spectrum, the model inherently exhibits CP violation with
$\delta$ approximately equal to $\pm 40^{o}$.

In summary, our investigation into the 2HDM model with GCP symmetries reduces
the number of free parameters, leading to precise predictions for neutrino
mixing angles and distinct CP-violating scenarios, shedding light on its
unique phenomenology.

Acknowledgments
Tapender acknowledges the financial support provided by Central University of
Himachal Pradesh. The authors, also, acknowledge Department of Physics and
Astronomical Science for providing necessary facility to carry out this work.

## References

  * [1] P. Minkowski, Phys. Lett. B 67, 421-428 (1977).
  * [2] T. Yanagida, Conf. Proc. C 7902131, 95-99 (1979).
  * [3] S. L. Glashow, NATO Sci. Ser. B 61, 687 (1980).
  * [4] M. Gell-Mann, P. Ramond and R. Slansky, Conf. Proc. C 790927, 315-321 (1979).
  * [5] R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980).
  * [6] T. D. Lee, Phys. Rev. D 8, 1226-1239 (1973).
  * [7] D. J. Castano, E. J. Piard and P. Ramond, Phys. Rev. D 49, 4882-4901 (1994).
  * [8] S. Weinberg, Phys. Rev. Lett. 40, 223-226 (1978).
  * [9] F. Wilczek, Phys. Rev. Lett. 40, 279-282 (1978).
  * [10] R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022).
  * [11] J. F. Gunion and H. E. Haber, Phys. Rev. D 67, 075019 (2003).
  * [12] S. L. Glashow and S. Weinberg, Phys. Rev. D 15, 1958 (1977).
  * [13] E. A. Paschos, Phys. Rev. D 15, 1966 (1977).
  * [14] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440-1443 (1977).
  * [15] S. Antusch, M. Drees, J. Kersten, M. Lindner and M. Ratz, Phys. Lett. B 525, 130-134 (2002).
  * [16] D. Atwood, S. Bar-Shalom and A. Soni, Phys. Lett. B 635, 112-117 (2006).
  * [17] J. D. Clarke, R. Foot and R. R. Volkas, Phys. Rev. D 92, no.3, 033006 (2015).
  * [18] Z. Liu and P. H. Gu, Nucl. Phys. B 915, 206-223 (2017).
  * [19] D. A. Camargo, M. D. Campos, T. B. de Melo and F. S. Queiroz, Phys. Lett. B 795, 319-326 (2019).
  * [20] D. Cogollo, R. D. Matheus, T. B. de Melo and F. S. Queiroz, Phys. Lett. B 797, 134813 (2019).
  * [21] L. Lopez Honorez, E. Nezri, J. F. Oliver and M. H. G. Tytgat, JCAP 02, 028 (2007).
  * [22] M. Gustafsson, E. Lundstrom, L. Bergstrom and J. Edsjo, Phys. Rev. Lett. 99, 041301 (2007).
  * [23] E. M. Dolle and S. Su, Phys. Rev. D 80, 055012 (2009).
  * [24] L. Lopez Honorez and C. E. Yaguna, JHEP 09, 046 (2010).
  * [25] G. Arcadi, Eur. Phys. J. C 78, no.10, 864 (2018).
  * [26] E. Bertuzzo, S. Jana, P. A. N. Machado and R. Zukanovich Funchal, Phys. Lett. B 791, 210-214 (2019).
  * [27] J. F. Gunion and H. E. Haber, Phys. Rev. D 72, 095002 (2005).
  * [28] P. M. Ferreira, H. E. Haber and J. P. Silva, Phys. Rev. D 79, 116004 (2009).
  * [29] P. M. Ferreira and J. P. Silva, Eur. Phys. J. C 69, 45-52 (2010).
  * [30] J. Velhinho, R. Santos and A. Barroso, Phys. Lett. B 322, 213-218 (1994).
  * [31] I. P. Ivanov, Phys. Rev. D 75, 035001 (2007) [erratum: Phys. Rev. D 76, 039902 (2007)].
  * [32] M. Maniatis, A. von Manteuffel, O. Nachtmann and F. Nagel, Eur. Phys. J. C 48, 805-823 (2006).
  * [33] C. C. Nishi, Phys. Rev. D 74, 036003 (2006) [erratum: Phys. Rev. D 76, 119901 (2007)].
  * [34] I. P. Ivanov, Phys. Rev. D 77, 015017 (2008).
  * [35] M. Maniatis, A. von Manteuffel and O. Nachtmann, Eur. Phys. J. C 57, 719-738 (2008).
  * [36] G. Ecker, W. Grimus and H. Neufeld, J. Phys. A 20, L807 (1987).
  * [37] P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola and J. W. F. Valle, JHEP 02, 071 (2021).
  * [38] C. Jarlskog, Phys. Rev. Lett. 55, 1039 (1985).
  * [39] S. M. Bilenky and S. T. Petcov, Rev. Mod. Phys. 59, 671 (1987) [erratum: Rev. Mod. Phys. 61, 169 (1989); erratum: Rev. Mod. Phys. 60, 575-575 (1988)].
  * [40] P. I. Krastev and S. T. Petcov, Phys. Lett. B 205, 84-92 (1988).
  * [41] A. Crivellin, A. Kokulu and C. Greub, Phys. Rev. D 87, no.9, 094031 (2013).


# Time-dependent density functional theory for X-ray near-edge spectroscopy

G.F. Bertsch  Institute for Nuclear Theory and Department of Physics,
University of Washington, Seattle, Washington  A. Lee  Department of Physics,
University of Washington, Seattle, Washington

###### Abstract

We derive a time-dependent density functional theory appropriate for
calculating the near-edge X-ray absorption spectrum in molecules and condensed
matter. The basic assumption is to increase the space of many-body wave
functions from one Slater determinant to two. The equations of motion derived
from Dirac’s variational principle provide an exact solution for the linear
response when the interaction Hamiltonian has only a core-electron field. The
equations can be solved numerically nearly as easily as the ordinary real-time
time-dependent Kohn-Sham equations. We carry out the solution under conditions
that permit comparison with the expected power-law behavior. Our extracted
power-law exponents are similar to those derived by Nozières and DeDominicis,
but are not in quantitative agreement. We argue that our calculational method
can be readily generalized to density functionals that take into account the
more general electron-electron interactions that are needed for treating
dynamic effects such as plasmon excitations.

###### pacs:

71.15.Mb,78.70.Dm,78.20.Bh

##  I Introduction

Time-dependent density functional theory (TDDFT) has proven to be a very
useful tool for calculating the linear response of condensed matter to
electromagnetic probes. The overall features of the dielectric function are
reproduced quite well, and the agreement at zero frequency in insulators is
often at the few percent level111It should be mentioned that there are also
well-known deficiencies in local density functionals, namely band gaps are too
high and excitons are missing from the theory.. It also describes the plasmon
peaks in the UV absorption spectrum and the corresponding energy loss spectrum
in inelastic electron scattering.

The good experience in TDDFT in the optical regime encourages its application
to X-ray absorption, which has a rich near-edge structure that is only
partially described by single-electron physics. Indeed, some aspects of the
X-ray response are accessible to TDDFTpr09 ; sa10 ; lo12 ; le12a . However,
other aspects, particularly the dynamics of core-hole relaxation, are beyond
the scope of the theory as presently formulated. This is because the
fundamental assumption that the wave function be represented by a single
Slater determinant is too restrictive when the core orbital is part of that
determinant. The single-determinant theory can describe correlations between
the core hole and the valence electron but not correlation effects within the
valence space. To overcome this deficiency, we propose here to extend the
TDDFT from one determinant to two. We derive the equations of motion in Sect.
II below. We refer the reader to Ref. ma12  for an review of the foundations
of TDDFT and its recent extensions. In Sect. III following we apply our
extension to the well-known Mahan-Nozières-De-Dominicis (MND) Hamiltonian. A
comprehensive discussion of analytic and numerical methods to solve is given
in the review by Ohtaka and Tanabe oh90 . The equations we solve numerically
are exact for the MND Hamiltonian. We will argue that they are much easier
than other methods to apply to DFT functionals which contain electron-electron
interactions.

##  II Extending TDDFT

We view the TDDFT as an approximation to Hamiltonian many-body theory taking
the wave function as a single Slater determinant (SD). The equations of motion
may be derived from the Dirac’s variational principle

$\delta\,\,\int dt\left<{\Psi(t)}\,\vrule\,{i{d\over
dt}-H}\,\vrule\,{\Psi(t)}\right>=0$ (1)

by varying the wave function $|\Psi\rangle$ in the space of Slater
determinants. The resulting time-dependent Kohn-Sham (KS) equations are then
solved for the time-dependent wave function $|\Psi(t)\rangle$. For discussion
of the action principle in the context of TDDFT, see Ref. vi08 ; ru12 .

The linear response to an operator $\cal O$ is obtained from the time-
dependent correlator $i\langle[{\cal O}(t),{\cal O}(0)]\rangle$ calculated
from

$R(t)=\lim_{\lambda\rightarrow 0}{1\over\lambda}\langle\Psi_{\lambda}(t)|{\cal
O}|\Psi_{\lambda}(t)\rangle.$ (2)

Here the initial state has been prepared by applying an impulsive field
$V(r,t)=\lambda{\cal O}\delta(t)$ to the KS ground state $|\Psi_{g}\rangle$.
In linear order, this modifies the initial wave function to

$|\Psi_{\lambda}(0_{+})\rangle=(1-i\lambda{\cal O})|\Psi_{g}\rangle.$ (3)

The resulting wave function is evolved by KS equations to determine the matrix
element in Eq. (2). The connection of Eq. (2) to the more familiar frequency-
dependent response $S(\omega)$ is given by

$S(\omega)\equiv\sum_{f}\langle f|{\cal
O}|0\rangle^{2}\delta(E_{f}-E_{0}-\omega)={1\over\pi}\int_{0}^{\infty}dtR(t)\sin\omega
t.$ (4)

Since $\lambda$ is small in Eq. (2), the evolved wave function is still
largely in the ground state with only a small amplitude of excited states.
This is fatal for calculating effects of the core hole excitation such as the
relaxation of the valence wave function in the presence of the core hole. By
considering separate determinants for the components of the wave function with
and without the core electron excitation, the correlations associated with
valence electrons can be treated as well as they are in the optical response.
There is no danger of violating the Fermi statistics because the two
components of the wave function are necessarily orthogonal. We note that
multicomponent TDDFT has already been derived as a generalization of the
Runge-Gross theorem and has been applied to a diatomic moleculele12b . The
multi-determinant time-dependent theory is also well-known in the literature
be00 ; al07 . However, it is typically based on a representation by particle-
hole excitations of a single determinant. Our derivation and equations of
motion are different, resembling more the multi-determinant theory of nuclear
excitations proposed in Ref. pu13 .

Our starting point is the following ansatz for the variational wave function,

$|\Psi\rangle=a_{g}c^{\dagger}_{c}|\Psi_{g}\rangle+a_{c}|\Psi_{c}\rangle=a_{g}c_{c}^{\dagger}\prod^{N_{e}}_{\alpha}c^{\dagger}_{\alpha}|\rangle+a_{c}\prod^{N_{e}+1}_{\beta}c_{\beta}^{\dagger}|\rangle,$
(5)

where $N_{e}$ is the number of active electrons in addition to the core
electron. Here $g,c$ in the middle equality designate the determinants
associated with the ground and core-excited wave functions, respectively. The
determinants are given more explicitly in the second equality, with
$c^{\dagger}_{\alpha},c^{\dagger}_{\beta}$ creation operators in the valence
band and $c^{\dagger}_{c}$ the creation operator for the core electron. The
two sets of valence-band orbitals are expressed in terms of the valence-band
basis states $i$ as $c^{\dagger}_{\alpha}=\sum a_{\alpha i}c^{\dagger}_{i}$
and $c^{\dagger}_{\beta}=\sum a_{\beta i}c^{\dagger}_{i}$. Typically, the
expansion coefficients $a_{\alpha i},a_{\beta i}$ satisfy the orthonormality
conditions $\sum_{i}a^{*}_{\alpha
i}a_{\alpha^{\prime}i}=\delta_{\alpha\alpha^{\prime}}$, etc. The amplitudes of
the two SD’s, $a_{g}$ and $a_{c}$, should satisfy the normalization condition
$|a_{g}|^{2}+|a_{c}|^{2}=1$.

The MND Hamiltonian has the form

$\hat{H}=\hat{H}_{v}(c^{\dagger}_{i},c_{i})+\hat{H}_{c}(c^{\dagger}_{i},c_{i})c_{c}c_{c}^{\dagger}.$
(6)

The first term is the valence Hamiltonian to be constructed from the
corresponding Kohn-Sham density functional. The second term adds the
excitation energy of the core hole as well as its field acting on the valence
electrons.

The variation in Eq. (1) is to be carried out with respect to changes in the
wave function $|\Psi\rangle$ that preserve its character as a sum of two SD’s.
In the single-determinant theory, one takes the variational derivatives of
$|\Psi\rangle$ with respect to $a_{\alpha,i}$ treating them as independent
variables. This results in the usual time-dependent Kohn-Sham (KS) with the
single-particle Hamiltonian given by222 The ordering of the operators in
$\hat{H}_{v}$ is responsible for the phase factor $s_{ij}=\pm 1$..

$\hat{H}_{KS}=\sum_{i,j}\langle\Psi|{\delta^{2}\hat{H}_{v}\over\delta
c^{\dagger}_{i}\delta c_{j}}|\Psi\rangle c^{\dagger}_{i}s_{ij}c_{j}$ (7)

However, as a consequence of the overcompleteness of the variables, the
overall phase of the wave function no longer has any physical meaning. For
example, if the orbitals are eigenstates of the KS Hamiltonian, the overall
phase is $\exp(-i\sum_{n}^{N_{e}}\epsilon_{n}t)$ where $\varepsilon_{n}$ are
the KS eigenvalues. The correct phase is $\exp(-i\langle\Psi|H|\Psi\rangle
t)$; the two are only equal in the absence of electron-electron interactions.
This phase plays no role in the single-determinant theory, but with two
determinants it is crucial to have correct relative phases.

The proper procedure to apply Eq. (1) is require that the wave function
variations $|\delta\Psi\rangle$ in derived equations

$\left<{\delta\Psi}\,\vrule\,{i{d\over
dt}-\hat{H}}\,\vrule\,{\Psi(t)}\right>=0.$ (8)

are independent of each other. An orthogonal (and thus independent) set of
wave functions may be defined by making use of Thouless’s representationRS  of
the SD’s. The equations of motion are obtained by taking $|\delta\Psi\rangle$
as the set of 1-particle 1-hole excitations of the instantaneous SD,

$|\delta\Psi\rangle\in|\alpha_{p}\alpha+h\rangle\equiv
c^{\dagger}_{\alpha_{p}}c_{\alpha_{h}}|\Psi(t)\rangle.$ (9)

in accordance with Thouless’s theorem. For a state of $N_{e}$ particles in a
basis of dimension $N_{b}$, there are $N_{e}(N_{b}-N_{e})$ particle-hole
amplitudes to be determined compared to the $N_{e}N_{b}$ amplitudes in the
representation Eq. (5). However, the use of Eq. (9) requires calculating both
particle and hole orbitals in an instantaneous basis, which is very costly in
carrying out the time evolution. An easier way to avoid the phase introduced
by the Kohn-Sham single-particle Hamiltonian is by projection. The action of
$\hat{H}_{KS}$ on the SD can be expressed in the instantaneous particle-hole
basis as

$\hat{H}_{KS}|\Psi\rangle=E_{KS}|\Psi\rangle+\sum_{\alpha_{p},\alpha_{h}}v(\alpha_{p},\alpha_{h})|\alpha_{p}\alpha_{h}\rangle$
(10)

where $E_{KS}=\langle\Psi|\hat{H}_{KS}|\Psi\rangle$. The unwanted first term
can be removed in any basis simply by updating the wave function using the
projected KS Hamiltonian $\hat{H}_{KS}-E_{KS}$. Thus, the single-particle
orbitals are calculated as usual, but the phase of the SD is corrected by
$\exp(+i\langle\hat{H}_{KS}\rangle\Delta t)$ at each time step. For our
numerical example below, the problem does not arise because there is no
electron-electron interaction in the valence space.

To summarize, we solve independently the time-dependent Kohn-Sham equations
for $|\Psi_{g}\rangle$ and $|\Psi_{c}\rangle$. The two determinants are
coupled by the X-ray photon interaction,

$\hat{H}_{x}=v_{x}(c^{\dagger}_{x}c_{c}+c^{\dagger}_{c}c_{x}),\,\,\,\,\,\,\,c^{\dagger}_{x}=\sum_{i}f_{i}c^{\dagger}_{i}$
(11)

with $f$ a form factor. Varying with respect to $a_{g},a_{c}$ one obtains the
$2\times 2$ matrix equation for these variables,

$i{d\over dt}\left(\begin{array}[]{c}a_{c}\\\
a_{g}\end{array}\right)=a_{c}\left(\begin{array}[]{cc}\langle\Psi_{c}|\hat{H}-E_{KS,c}|\Psi_{c}\rangle&v_{x}\langle\Psi_{c}|\hat{H}_{x}c^{\dagger}_{c}|\Psi_{g}\rangle\\\
v_{x}\langle\Psi_{g}|c_{c}\hat{H}_{x}|\Psi_{c}\rangle&\langle\Psi_{g}|\hat{H}-E_{KS,g}|\Psi_{g}\rangle\\\
\end{array}\right)\left(\begin{array}[]{c}a_{c}\\\ a_{g}\end{array}\right).$
(12)

The hermiticity of $\hat{H}_{x}$ ensures that the normalization condition
remains satisfied during the course of the evolution. The off-diagonal matrix
element in this equation is expressible as the $N\\!+\\!1\times N\\!+\\!1$
determinant

$\langle\Psi_{c}|c^{\dagger}_{x}|\Psi_{g}\rangle=\left|\begin{array}[]{ccc}\langle\beta_{1}|\alpha_{1}\rangle&...&\langle\beta_{1}|x\rangle\\\
\langle\beta_{2}|\alpha_{1}\rangle&...&\\\ ...&&\\\
\langle\beta_{N+1}|\alpha_{1}\rangle&...&\langle\beta_{N+1}|x\rangle\\\
\end{array}\right|.$ (13)

While this determinant is well-known in the analytic theory of the near-edge
responseoh90 , it is absent from the usual time-dependent Kohn-Sham theory
based on a single Slater determinant.

To evaluate the linear response to the field of the X-ray photon, we start
with the ground state wave function at time zero, $|\Psi_{g}\rangle$. We now
perturb the system by an impulsive X-ray field, $\lambda\hat{H}_{x}\delta(t)$.
The immediate evolution introduces a small component of the core-excited state
into the wave function,

$|\Psi(0_{+})\rangle=c^{\dagger}_{c}|\Psi_{g}(0)\rangle+i\lambda
v_{x}|\Psi_{c}(0)\rangle,$ (14)

where $|\Psi_{c}(0)\rangle=c^{\dagger}_{x}|\Psi_{g}(0)\rangle$. Eq. (14) has
the required form as a sum of two determinants. Each is evolved in time with
its own Kohn-Sham Hamiltonian. Then the real-time response from Eq. (2) is

$R(t)=2v_{x}^{2}{\rm
Re}~{}\langle\Psi_{c}(t)|c^{\dagger}_{x}|\Psi_{g}(t)\rangle.$ (15)

This can be easily Fourier-transformed by Eq. (4 ) to give the absorption
spectrum.

Our procedure provides an exact solution for the linear response if
$\hat{H}_{v}$ and $\hat{H}_{c}$ are strictly one-body operators. This is
because the intrinsically two-body part of the Hamiltonian does not make two-
particle excitations or entangle the two Slater determinants after the
initialization.

##  III Numerical calculations

In this section we demonstrate the practicality of the method as applied to
the MND Hamiltonian. The computer codes employed here are available at Ref.
computer .

We write the two terms in the Hamiltonian as

$\hat{H}_{v}={E_{b}\over
N_{b}-1}\sum_{i=1}^{N_{b}}(i-N_{b}/2)c^{\dagger}_{i}c_{i}$ (16)

and

$\hat{H}_{c}={v_{c}\over
N_{b}}c^{\dagger}_{x}c_{x},\,\,\,\,c^{\dagger}_{x}=\sum_{i}^{N_{b}}c_{i}^{\dagger}.$
(17)

Here $E_{b}$ is the width of the band, and $N_{b}$ is thenumber of orbitals in
the band. We shall express energies in units of $E_{b}$, and time in units of
$\hbar/E_{b}$. We start with a half-filled band, taking the number of valence
electrons $N_{e}$ to be $N_{e}=N_{b}/2$. We present calculations for the
parameter sets listed in Table I.

As explained in the literature oh90 , the core-hole interaction strength
$v_{c}$ is not the most physically direct quantity determining the near-edge
response. The effect of the Fermi-surface edge is more closely related to the
shift of the single-particle orbital energies due to the core hole. Calling
the shift $\Delta\varepsilon$, the relevant quantity is

$\Delta\varepsilon{dn\over d\varepsilon}={\delta\over\pi}$ (18)

where $dn/d\varepsilon$ is the density of orbital states at the Fermi level.
In the last equality, this is related to the scattering phase shift $\delta$
at the Fermi surface. The values of $\delta$ associated with the computed
parameter sets are given in the last column of Table I.

The Green’s function theory in Ref. no69  for the time-dependent response
decomposes it into two factors, the overlap of Fermi sea determinants and the
Green’s function of the electron that was promoted to the valence band. We
write the overlap of the Fermi sea determinants as

$G(t)=\langle\Psi_{g}|e^{-i(\hat{H}_{v}+\hat{H}_{c})t}|\Psi_{g}\rangle.$ (19)

The main quantity of interest is the determinant in Eq. (15) which we call
$g_{c}$, as

$g_{c}(t)=\langle\Psi_{c}(t)|c^{\dagger}_{x}\Psi_{g}(0)\rangle$ (20)

Nozières and De Dominicis decompose it into two factors, $g_{c}(t)=g(t)G(t)$.
We will not make use of that separation.

###  III.1 Fermi sea evolution

We first examine the Fermi sea overlap. To remove the phase of the core-
excited ground state, we will examine the quantity

$G^{\prime}(t)=e^{i\sum_{\alpha}\varepsilon_{\alpha}t}\langle\Psi_{g}|e^{-i(\hat{H}_{v}+\hat{H}_{c})t}|\Psi_{g}\rangle,$
(21)

where $\epsilon_{\alpha}$ are the Kohn-Sham eigenvalues of the ground-state
orbitals. Figure 1 shows ${\rm Re}~{}G^{\prime}(t)$ for parameter set A of
Table I.

Figure 1:  ${\rm Re}~{}G^{\prime}(t)$ as a function of time for parameter set
A. The line shows a visual power-law fit, $G^{\prime}(t)\sim t^{-0.13}$.

It is plotted on a log-log scale to facilitate comparison with result of Ref.
no69 ,

$G(t)\sim(1/t)^{(\delta/\pi)^{2}}.$ (22)

The predicted power-law dependence should be applicable over the time domain
starting from $t_{0}\sim 1/E_{b}$ and going to $t_{1}\sim
dn/d\varepsilon=N_{b}/E_{b}$, the time necessary to resolve individual
orbitals in the band. In our units the range is $(t_{0},t_{1})=(1,N_{b})$. One
notices immediately that $G^{\prime}(t)$ has a considerable oscillatory
component. The oscillation has been found in other treatments of the problem
as well, eg. (oh86, , Eq. (3.4)). As discussed in Ref. oh90 , the oscillation
may be attributed to the deeply bound orbital at the bottom edge of the
valence band.

case | $N_{b}$ | $N_{e}$ | $v_{c}$ | $\delta/\pi$
---|---|---|---|---
A | 256 | 128 | -0.8 | 0.38
B | 8 | 4 | -0.8 | 0.39
C | 512 | 256 | -0.8 | 0.38
Z | 256 | 128 | 0 | 0
Table 1: Parameter sets for the Hamiltonian Eq. (16,17).

The line in the graph corresponds to the power law $G(t)\sim(1/t)^{\gamma}$
with $\gamma=0.13$. This is rather close to the predicted power law derived in
Ref. no69 , $\gamma=(\delta/\pi)^{2}\approx 0.14$.

The spectral function associated with $G^{\prime}(t)$ is its Fourier
transform,

$G^{\prime}(\omega)=\int_{0}^{\infty}dte^{i\omega t}G^{\prime}(t).$ (23)

Fig. 2 plots ${\rm Re}~{}G^{\prime}(\omega)$ for parameter sets A and B. For
the set A shown in the left-hand panel one can see the peaks associated with
individual states of the many-particle wave function. The dimension of the
many-particle space is given by ${N_{b}\choose{N_{e}}}=20$. The ground state
is the peak on the far left, and 8-9 other states are visible in the plot. The
right-hand plot shows ${\rm Re}~{}G^{\prime}(\omega)$ for parameter set A.
Here the individual states are so closely spaced that one can see only smooth
curves. There are clearly two peaks in the spectrum, one associated with the
ground state and its low-energy excitations, and the other associated with a
localized orbital bound or nearly bound orbital to the core hole. In Fig. 3 we
have replotted the ${\rm Re}~{}G^{\prime}(\omega)$ ground-state peak on a log-
log scale to make visible a power-law dependence on $\omega$. The expected
range of validity for a power law is within the interval
$(\omega_{0},\omega_{1})=(1/N_{b},1)$ in our units. The line in Fig. 3 shows
the the power-law $\omega^{-1.13}$. We can see that it provides a reasonable
fit in the range $(0.02,0.3)$ with some oscillation at low frequency.

Figure 2:  ${\rm Re}~{}G^{\prime}(\omega)$ as a function of $\omega$. Left-
hand panel: results for a low-dimensional system, parameter set B. Right-hand
panel: results for parameter set A.  Figure 3:  The same $G^{\prime}(\omega)$
as in Fig.2b, plotted on a log-log scale. The line shows a visual power-law
fit, $G^{\prime}(\omega)\sim\omega^{-1.13}$.

###  III.2 Inclusion of the core electron

We now examine the propagation of the core-hole excited determinant with the
core electron promoted to the valence band. The number of electrons in the
determinant is now $N_{b}/2+1$. The initial wave function has equal amplitudes
for the x electron in all the unoccupied orbitals; it thus has the same
localization as the core-hole potential. Just as a reminder of the non-
interacting physics, we show in Fig. (4) the imaginary part of $g_{c}(\omega)$
at $v_{c}=0$. It is uniform across the region of unoccupied orbitals, with
sharp edges at the Fermi surface and at the top of the band.

Figure 4:  Left panel: core-excited Green’s function with no interaction,
parameter set Z. Right panel: the same quantity with parameter set A. Note the
difference in vertical scale.

The right-hand panel shows the Green’s function with parameter set A. Note
that the peak associated with a hole at the bottom of the valence band is
missing. Evidently, the electron added by $c^{\dagger}_{x}$ operator ensures
with a high probability that the hole is filled. The results in the right-hand
panel are plotted in Fig. 5 on a log-log scale. The line is a visual fit to
power-law behavior. Evidently, a power law gives a reasonable description over
the energy interval $0.02-0.3$. According to the theory, the exponent is
determined by the phase shift according to the dependence 333This is for a
single partial wave and spin projection..

$\gamma=-2{\delta\over\pi}+{\delta^{2}\over\pi^{2}}.$ (24)

Taking $\delta/\pi$ from Table I, the predicted value is $\gamma=-0.62$,
compared to $\gamma=-0.85$ from the fit.

The small disagreement we find here persists over a large range of $\delta$.
Figure 6. shows a comparison over the range of $\delta$ accessible to the
Hamiltonian. We see that the exponent is proportional to $\delta$ for small
$\delta$ as in Eq. (24), but the coefficient is somewhat higher.

Figure 5:  The determinant $g_{c}(\omega)$, plotted on a log-log scale. The
line shows an approximate power law fit, $g_{c}(\omega)\sim\omega^{-0.85}$.
Figure 6:  Solid circles: the power-law exponent extracted from the extended
TDDFT response for the space $(N_{b},N_{e})=(256,128)$. Dashed line: the
analytic formula Eq. (24). The numerically computed exponent was extracted
from the calculated $g_{c}(\omega)$ at $\omega=0.03$ and $0.2$ .

###  III.3 Other numerical treatments

Various numerical treatments of the MND Hamiltonian have been discussed and
reviewed in Ref. oh90 . The two main approaches are the Green’s function
formulationsc77 ; gr77 ; pr01 ; ba80 ; oh86  and the formulation in terms of
many-body determinantal wave functionsko73 ; sw79 ; do80 ; fe80 ; ma80 . The
former requires constructing functions of at least two variables in the time
or frequency domain, governed by equations that are nonlocal in those
variables. In this respect, the real-time wave function method is much more
efficient since there is only one time variable and the equations to be solved
are local in time. An early numerical work following the wave function
approach is Ref. ko73 . These authors constructed the numerically exact
solutions of the core-hole excited Hamiltonian, and then used the ground state
and one-particle excitations of the Hamiltonian as the final states. This
procedure of enumerating the eigenstates was also used in Refs. fe80 . The
wave function approach was also used in Ref. ma80 , and the determinant was
evaluated in real time, as in our approach. However, the procedure adopted
there was based on the formulation of Ref. co71  which requires a matrix
inversion. The near-singularity of the matrix apparently caused numerical
difficulties that do not arise in the real-time approach. From a computational
point of view, our approach is closest to that used in Ref. sw79  and do80 .
We note that these authors found that the critical exponents of the analytic
Green’s function treatment were only in qualitative agreement with the
numerical results outside of a very small interval near $\omega=0$.

##  IV Summary and outlook

We have derived an extension of time-dependent density functional theory that
contains at least some of the subtle many-particle physics of X-ray near-edge
absorption in metals. Numerically, the real-time theory is easy to carry out
if the time-dependent electron-electron interaction is neglected. The
absorption shows that the X-ray absorption power-law behavior is in
qualitative agreement with the analytic results of Ref. ro69 , but not
identical to them. A similar conclusion was obtained in Ref. sw79 .

Physically, the most important effect of the many-electron physics is core-
hole screening. There are several numerical calculations in the literature
that follow the Green’s formalism of Ref. no69  and focus on this screening
effect. For the absorption spectra, a commonly used approximation treats the
system as fully relaxed in the presence of the core hole. Good fits can also
be obtained under the assumption that the dynamic screening reduces the core-
hole effects by a factor of twole12a . That work also presented a real-time
dynamic calculation, but apparently used a diagonal approximation to Eq. (15).
In any case, dynamic effects related to the core hole can be easily calculated
in the real time method, so there is no reason to use any of these
approximations.

So far we have not discussed the electron-electron interactions within the
valence band. They are potentially important and are needed to treat the
additional screening associated with the plasmon degree of freedom. Langreth
has proposed a way to include plasmon effects in the Green’s approach la70
and it was applied with some success in Ref. gu12 . However, it involves
distinct calculations for the plasmon physics and for the X-ray absorption. In
contrast, the real-time TDDFT provides a unified framework for including the
electron-electron interaction in the calculation. In the two-determinant
theory one can simply add the Coulomb field of the instantaneous charge
density of $|\Psi_{c}\rangle$ to the field generated by $v_{c}$. Of course,
the presence of the interaction requires considerably more computational
effort. The Coulomb field has to be calculated at each time step. Also, the
overall phase of the $|\Psi_{c}\rangle$ has to be computed using one of the
methods discussed in Sect. II. We believe that the problem is still
computationally quite tractable; we leave the implementation to a future
publication.

## Acknowledgment

We acknowledge discussions with J.J. Rehr, A. Bulgac, and K. Yabana. We also
thank U. Mosel for bringing Ref. pu13  to our attention. This work was
supported by the National Science Foundation under Grant PHY-0835543 and by
the (US) Department of Energy under Grant No. DE-FG02-00ER41132.

## References

  * (1) M. Prange, et al., Phys. Rev. B 80 155110 (2009).
  * (2) A. Sakko, A. Rubio, M. Hakala, and K. Hamalainen, J. Chem. Phys. 133 174111 (2010).
  * (3) K. Lopata, B. Van Kuiken, M. Khalil, and N. Govind, J. Chem. Theory Comput. 8(9) 3284 (2012).
  * (4) A. Lee, F. Vila, and J. Rehr, Phys. Rev. B 86 115107 (2012).
  * (5) M. Marque, N. Maitra, F. Nogueira, and E. Gross, eds. “Fundamentals of time-dependent density-functional theory,” Lecture Notes in Physics, 837 (2012).
  * (6) K. Ohtaka and Y. Tanabe, Rev. Mod. Phys. 62 929 (1990).
  * (7) G. Vignale, Phys. Rev. A 77 062511 (2008).
  * (8) M. Ruggenthaler and R. van Leeuwen, Ref. ma12 , pp. 206-210.
  * (9) R. van Leeuwen and E. Gross, Ref. ma12 , p. 249.
  * (10) M. Beck, et al., Phys. Rep. 324 1 (2000)
  * (11) O. Alon, J. Chem. Phys. 127 154103 (2007).
  * (12) G. Puddu, arXiv:1208.0122v3, (2013).
  * (13) P. Ring and P. Schuck, “The Nuclear Many-Body Problem,” (Springer, 1980).
  * (14) P. Nozières and C.T. De Dominicis, Phys. Rev. 178 1097 (1969).
  * (15) T. Privalov, R. Gel’mukhanov, and H. Agren, Phys. Rev. B64 165115 (2001).
  * (16) K. Ohtaka and Y. Tanabe, Phys. Rev. B 34 3717 (1986).
  * (17) K. Schn̈hammer and O. Gunnarsson, Solid State Comm. 23 691 (1977).
  * (18) V. Grebennikov, Y. Babanov, and O. Sokolov, Phys. Status Solidi B 80 73 (1977).
  * (19) U. von Barth and G. Grossmann, Physica Scripta 21 580 (1980).
  * (20) A. Kotani and Y. Toyozawa, J. Phys. Soc. Jpn. 35 1082 (1973).
  * (21) C. Swartz, J. Dow, and C. Flynn, Phys. Rev. Lett. 43 158 (1979).
  * (22) J.Dow and C. Flynn, J. Phys. C 13 1341 (1980).
  * (23) L. Feldkamp and L. Davis, Phys. Rev. B 22 4994 (1980).
  * (24) G. Mahan, Phys. Rev. B 21 1421 (1980).
  * (25) M. Combescot and P. Nozières, J. Phys. (Paris) 32 913 (1971).
  * (26) B. Roulet, J. Gavoret, and P. Nozières, Phys. Rev. 178 1072 (1969).
  * (27) D.C. Langreth, Phys. Rev. B1 471 (1970).
  * (28) M. Guzzo, et al., Eur. Phys. J. B 85 324 (2012).
  * (29) See the URL http://www.phys.washington.edu/users/bertsch/rt-xfs.tar or “Supplementary Material” accompanying the published journal article.