# K-essence scalar field as dynamical dark energy

L. A. García<EMAIL_ADDRESS>J. M. Tejeiro<EMAIL_ADDRESS>L.
Castañeda<EMAIL_ADDRESS>Observatorio Astronómico Nacional,
Universidad Nacional de Colombia.

###### Abstract

We study an early dark energy (EDE) model as a K-essence scalar field in the
framework of FLRW universe using an effective parametrization of the state
equation as a function of the redshift $z$ with the tracker condition during
radiation domination, but also demanding an accelerated expansion of the
universe at late times emulating cosmological constant. We found all the
dynamical variables of the EDE system. We use the luminosity distances of the
SNIA to get the best estimations for the free parameters of the model and
also, we constrain the model using primordial abundances of light nuclei in
BBN theory. We summarize the necessary conditions to achieve BBN predictions
and the accelerated expansion of the universe at late times.

###### pacs:

26.35.+c, 95.36.+d.

##  I Introduction

Recent observations of the luminosity distances of the SNIa perlmutter  reveal
that the expansion of the universe is accelerated and there is an unknown
matter–energy contribution as about 70% of the critical density, which is
smooth and has negative pressure. In order to explain this phenomenon, it has
been proposed many plausible solutions: the cosmological constant $\Lambda$
which is related to the vacuum energy of the quantum fields carroll ,
Quintessence fields (with the state equation
$\omega=\frac{p_{Q}}{\rho_{Q}}$=constant), Kessence, Taquionic fields,
frustrated topological defects, extra–dimensions, massive (or massless)
fermionic fields, galileons, effective parametrizations of the state equation,
primordial magnetic fields, holographic models, etc. All these proposals have
been used to model this contribution predicted by the Friedmann equations in
the framework of the General Relativity. Also, there are other possible
options as Modified Gravity, where the accelerated expansion effect is
geometric and not as a matter–energy form.
Nowadays, the current paradigm is the $\Lambda$CDM model, which is so far the
best fitting to the present observations, even though the conceptual problems
that persist with the nature of $\Lambda$.

On the other hand, primordial abundances of the light nuclei that were formed
during Big Bang Nucleosynthesis (BBN) are well observed and quantified with
astrophysical methods, specially the mass fraction of the ${}^{4}He$.
Actually, there is an extended theory of BBN, introduced by Alpher, Bethe and
Gamov alpher ; gamov , and it has been developed many numerical codes that
resolve the Boltzmann equation for each isotope including the cosmological
background (based on Kawano code). There is a good agreement between the
predicted abundances and the observations, however, the main problem remains
with ${}^{7}Li$. Despite the efforts that have been made, including
corrections on the cross sections of this element, the discrepancies are not
negligible.

With the aim to enhance the primordial abundances calculated from BBN and
also, give a plausible explanation of the accelerated expansion of the
universe, we have proposed a model of dark energy which has a non–null
contribution at early times to increase the Hubble radius during radiation
domination era and influence the Boltzmann equations that determine the
evolution of the light abundances. All the conditions that allow us to
describe the early dark energy are achieved with a K-essence scalar field,
which is characterized by its state equation that overcomes the attractors
defined by a dynamical system.

In the first section, we expose the main conditions imposed to the scalar
field and we resolve the dynamical system that appears from the cosmological
assumptions. At the second part, we propose the effective parametrization of
the state equation, determinate the best estimations of the free parameters of
the model using the luminosity distances of the SNIa from the Union dataset
and derive analytical expressions for the dynamical variables of the K-essence
system. In the IV section, we test our model with some standard proves and
summarize the main differences between this model and $\Lambda$CDM.

Taking into account that there is a non-negligible energy density of the
K-essence field during radiation era, we compute the BBN abundances, including
the scalar field degrees of freedom in Hubble parameter, following the
derivation and discussion made by Bernstein bernstein2  and using the
available codes for Nucleosynthesis. Finally, we find the values of
$\Omega_{B}$ and $\eta_{B}$ predicted by our model as a result of the BBN
calculations.

##  II K-essence scalar field and the dynamical system

The K-essence scalar fields appeared at the late 90s with the K-inflation
model proposed by Armendariz-Picon armendariz1 ; armendariz2 . However, the
idea was extended to describe a dynamical dark energy contribution, taking
into account that this field can track during radiation domination epoch and
also, it could avoid the fine–tuning of the initial values of the field and
its velocity. These features are well known for the system, so the challenge
is resolve the evolution equations of the field in the FLRW spatially flat
universe and without cosmological constant.

The lagrangian of the K-essence scalar field is given by:

$p(X,\phi)=K(\phi)L(X),$ (1)

where $X=-\frac{1}{2}\nabla^{\alpha}\phi\nabla_{\alpha}\phi$ is the kinetic
energy of the field and $v$ its velocity $v=\frac{d\phi}{dt}=\sqrt{-2X}>0$.
Rewriting (1) in terms of $v$, the lagrangian has the following form :

$p(v,\phi)=K(\phi)Q(v).$ (2)

The most general action that describes the K-essence field in a cosmological
plasm is given by:

$S=\int d^{4}x\sqrt{-g}\left(\frac{R}{2\kappa^{2}}+p(\phi,X)\right)+S_{B},$
(3)

with $\kappa^{2}=8\pi G$ and $S_{B}$ is the action of the background matter.
The signature that is used is $\\{-1,+1,+1,+1\\}$ The equations of motion of
the field come from the variation of the lagragian with respect to the field
$\phi$:

$\left(\frac{1}{v}\frac{\partial p}{\partial v}+v\frac{\partial}{\partial
v}\left(\frac{1}{v}\frac{\partial p}{\partial
v}\right)\right)\ddot{\phi}+\frac{\partial p}{\partial
v}(3H)+v\frac{\partial^{2}p}{\partial\phi\partial v}-\frac{\partial
p}{\partial\phi}=0.$ (4)

From (3), it is obtained the energy–momentum tensor:

$T_{\mu\nu}=\partial_{\mu}\phi\partial_{\nu}\phi-p(\phi,v)g_{\mu\nu}.$ (5)

Since the scalar field can be described as a perfect fluid, the energy density
and the pressure are defined by:

$\rho_{\phi}=K(\phi)\left(v\frac{\partial Q}{\partial v}-Q\right),$ (6)

$p_{\phi}=K(\phi)Q(v).$ (7)

In addition, the adiabatic velocity of sound for the K-essence field is given
by:

$C_{s}^{2}=\dfrac{Q^{\prime}}{vQ^{\prime\prime}}.$ (8)

where $\prime$ denotes derivative with respect to the velocity of the field
$v$. $C_{s}^{2}$ gives relevant information of the stability of the
perturbations associated with the K-essence field.
In order to resolve the (4) and find an explicit form of (6) and (7), there
have been suggested many alternatives: fixing a specific function of $\phi(t)$
or $v(t)$ chiba , making redefinition of the field to face a modified
lagrangian $\bar{Q}(v)$ copeland , considering pure kinetic K-essence model
stiff ; linder ; nojiri ; yang ; scherrer2  or imposing Slow roll conditions
on $K(\phi)$.

However, we want to resolve the complete dynamical system defined by the
Friedmann and the continuity equations for non-interactive fluids in the
cosmological background: a matter (or radiation component) and the K-essence
scalar field.

$H=\frac{\dot{a}}{a}=\kappa\sqrt{\rho_{\phi}+\rho_{m}},\>\>\>\>\>\>$ (9)
$\frac{\ddot{a}}{a}=-\frac{\kappa^{2}}{6}((1+3\omega_{\phi})\rho_{\phi}+(1+3\omega_{m})\rho_{m}),\>\>\>\>\>\>$
(10)
$\dot{\rho_{m}}=-3H(1+\omega_{m})\rho_{m},\>\>\>\>\>\>\>\>\>\>\>\>\dot{\rho_{\phi}}=-3H(1+\omega_{\phi})\rho_{\phi},\>\>\>\>\>\>$
(11)

Reexpressing the equiations (9), (10), (11) in terms of the velocity and the
dimensionless variable $F=\frac{\rho_{m}}{\rho_{\phi}+\rho_{m}}$ (matter
energy density fraction):

$\frac{dv}{d\phi}=-C_{s}^{2}\left[\frac{(lnK)_{,\phi}v}{1+\omega_{\phi}}+3\kappa\left(\frac{K\tilde{\rho_{\phi}}}{1-F}\right)^{1/2}\right],$
(12)

$\frac{dF}{d\phi}=-\frac{3\kappa}{v}F\sqrt{1-F}\sqrt{K\tilde{\rho_{\phi}}}\left(\omega_{m}-\omega_{\phi}\right).$
(13)

In addition, it must be fulfilled the following conditions:

$\frac{d\tilde{\rho_{\phi}}}{dv}=\frac{(1+\omega_{\phi})}{vC_{s}^{2}}\tilde{\rho_{\phi}},$
(14)

$\frac{d\omega_{\phi}}{dv}=\frac{1+\omega_{\phi}}{v}\left(1-\frac{\omega_{\phi}}{C_{s}^{2}}\right).$
(15)

It is assumed an asymptotic behaviour for the function $K(\phi)$ kang :

$K(\phi)=\frac{1+K_{0}(\phi)}{\phi^{2}},\>\>\>\>\>\>\>\>\>\lim_{\phi\to\infty}K_{0}(\phi)=0.$
(16)

To reach an attractor during radiation domination epoch, the dynamical system
$\\{v(\phi),R(\phi)\\}$ must fulfill the ansatz:

$v(\phi)=v_{0}-A(\phi),\>\>\>\>\>\>\>\>\>F(\phi)=F_{0}-B(\phi).$ (17)

where $A(\phi),B(\phi)\rightarrow 0$ monotonically for $\phi\to\infty$ (or
equivalently, $v=v_{rad}$ and $F\sim 0$). The ansatz implies the following
physical conditions on the system:

$\omega_{m}=\omega_{\phi}(v_{0}),\>\>\>\>\>\>\>\>\tilde{\rho_{\phi}}(\phi)\neq
0,\>\>\>\>\>\>\>\>C_{s}^{2}(v_{0})>\omega_{m}$ (18)
$\omega_{\phi}^{\prime}(v_{0})=\left(1-\frac{\omega_{m}}{C_{s}^{2}(v_{0})}\right)\frac{1+\omega_{m}}{v_{0}}\neq
0.$ (19)

There is also a De-Sitter attractor given by the condition $R\sim 0$, which
guarantees the existence of a accelerated expansion at late times:

$v(\phi)=v_{s}-A(\phi),\>\>\>\>\>\>\>\>\>R(\phi)=B(\phi).$ (20)

The last ansatz entails a condition on the state equation on the vicinity of
the De-Sitter attractor:

$\dfrac{{1-\omega_{\phi}(v_{s})}}{{1+\omega_{\phi}(v_{s})}}>0,\>\>\>\>\>\>\>\>|\omega_{\phi}(v_{s})|<1.$
(21)

##  III Effective parametrization of $\omega_{\phi}$

In order to get a general solution for the dynamical system defined by (14),
(15), (18), (19) and (21), it has been proposed an effective parametrization
of the state equation from $z<10^{15}$ given by:

$\omega_{\phi}(z)=\frac{4/3}{\left(\frac{1+z_{d}}{1+z}\right)^{m}+1}-1,$ (22)

where $m$ is factor that modules the transitions between the attractors,
$z_{d}$ is a redshift in matter domination epoch defined by
$z_{d}=\dfrac{z_{eq}+z_{*}}{2}$ and $z_{*}$, the redshift where the De-Sitter
domination -accelerated expansion- begins.

The parametrization (22) respects all the conditions previously mentioned,
hence it is possible to resolve the functions related to the K-essence
lagrangian.
The energy density of the field $\rho_{\phi}$ (6) results:

$\int_{\rho}^{\rho_{0}}\frac{d\rho^{\prime}}{\rho^{\prime}}=-3\int_{a}^{1}\frac{(1+\omega_{\phi}(a^{\prime}))}{a^{\prime}}da^{\prime},$
(23)

integrating (23), it is obtained:

$\rho=\rho_{0}\cdot(1+z)^{4}\left[\frac{\left(\left(\frac{1+z_{d}}{1+z}\right)^{m}+1\right)}{\left(\left(1+z_{d}\right)^{m}+1\right)}\right]^{4/m}=\rho_{0}\cdot
f(z).$ (24)

with

$\displaystyle f(a)$
$\displaystyle=exp\left[-3\int_{a}^{1}\frac{(1+\omega_{\phi}(a^{\prime}))}{a^{\prime}}da^{\prime}\right]$
(25)
$\displaystyle=a^{-4}\left[\frac{\left(\left(\frac{a}{a_{d}}\right)^{m}+1\right)}{\left(\left(\frac{1}{a_{d}}\right)^{m}+1\right)}\right]^{4/m},$
(26)

Meanwhile, the fraction of the dark energy density
$\Omega_{\phi}=1-F=\frac{\rho_{\phi}}{\rho_{cr}}$:

$\Omega_{\phi}=\frac{\Omega_{\phi 0}\cdot f(a)}{\Omega_{\phi 0}\cdot
f(a)+\Omega_{m0}\cdot a^{-3}},$ (27)

The formal solution of (22) will be obtained with the best estimation of the
free parameters of the model $\\{\Omega_{\phi_{0}},m,z_{*}\\}$. For this
reason, we use the luminosity distances of the SNIa from the survey Supernova
Cosmology Project with $z>0.8$ to minimize the function $\chi^{2}$:

$\chi^{2}=\sum_{i=1}^{N\sim 59}\frac{[\mu_{i}-\mu(z_{i})]^{2}}{\sigma^{2}}$
(28)

with the distance modulus given by the expression
$\mu=m-M=5(log_{10}d_{L}(z)-1)$ and the luminosity distance of our model:

$d_{L}(z)=\frac{c(1+z)}{H_{0}}\int_{0}^{z}\frac{dz^{\prime}}{B(z)}.$ (29)

$B(z)=(\Omega_{\phi_{0}}f(z^{\prime};m,z_{*})+(1-\Omega_{\phi_{0}})(1+z^{\prime})^{3})^{1/2}$

(28) must be resolved together with the constraints:

$\sum_{i=1}^{N\sim
59}\left(\frac{\mu_{i}-\mu(z_{i})}{\sigma^{2}}\right)\left(\frac{\partial\mu(z_{i};\Omega_{\phi_{0}},m,z_{*})}{\partial
x}\right)=0,$ (30)

where $x=\Omega_{\phi_{0}},m,z_{*}$. In addition, we have imposed other 2
conditions for the parameters: the deceleration parameter has to be zero at
$z_{*}$, therefore $q(z_{*})=0$, but also $f(z\sim z_{BBN}\approx 10^{9})$ has
to overcome the maximum value at the Primordial Nucleosynthesis to contribute
with some relativistic degrees of freedom and enhance the predicted primordial
abundances.

Resolving (28) simultaneously as (30), it is found the following values for
the free parameters:

| K-essence (22)
---|---
$\Omega_{\phi_{0}}$ | 0.69
$m$ | 1.0
$z_{*}$ | 1.48
$\omega_{0}$ | \- 0.99

In the figure 1 is plotted the evolution of the state equation $\omega$ as
function of $a$. The field emulates radiation during this epoch and then it
evolves to the next attractor: De-Sitter. Here is clear the tracker behaviour
imposed in the dynamical system. In addition, it is shown that the field is
relaxing to the asymptotic $\Lambda$CDM model for late times (during De-Sitter
attractor.)

Figure 1: State equation of K-essence field as a function of $a$.

Figure 2 shows the behaviour of the the function $f(z)$ that characterizes the
dark energy density evolution in the model. During radiation domination epoch,
the field scales as radiation $\rho\propto(1+z)^{4}$ until $z_{eq}$. After
that, $\rho$ has a complex behaviour which guarantees the second attractor
will be reached. At this point, the K-essence scalar field evolves in the De-
Sitter attractor and its state equation goes asymptotically to $-1$ (as
cosmological constant).

Figure 2: Function $f(z)$: dark energy density evolution.

The figure 3 displays the luminosity distance for the model compared with the
predicted by $\Lambda$CDM. The shift between the curves shows that the
luminosity distance is upper than the associated with the $\Lambda$CDM model,
because the matter density today predicted by our model is higher than the
second one (compared with WMAP-7 $\\{\Omega_{\phi
0},\Omega_{m0}\\}=\\{0.734\pm 0.029,0.266\pm 0.029\\}$
111lambda.gsfc.nasa.gov/product/map/dr4/parameters.cfm.

However, for low redshifts the luminosity distance grows linearly independent
on the model at this regime.

Figure 3: Luminosity distance as a function of $z$.

On the other hand, figure 4 displays the evolution of the dark energy density
fraction for the model with the estimations of the free parameters (III):

Figure 4: Evolution of the dark energy density fraction as a function of $z$.

Up to now, we are interested to resolve the dynamical variables of the
K-essence according to the effective parametrization we have proposed. The
first one of these quantities is the adiabatic velocity of the sound
$C_{s}^{2}$, which has a behaviour defined by equation (31).

$C_{s}^{2}=\omega_{\phi}+\frac{m}{4}\left(\frac{1+z_{d}}{1+z}\right)^{m}(\omega_{\phi}+1).$
(31)

Using (22), we complete the evolution of $C_{s}^{2}$. The behaviour of the
adiabatic velocity is plotted in the figure 5.

Figure 5: Adiabatic sound velocity evolution within the model.

$C_{s}^{2}$ fulfills the condition (18) in the radiation domination epoch
attractor. Whatsoever, there is not a straightforward physical interpretation
of the adiabatic velocity during De-Sitter attractor, because this value
implies a complex Young module of the plasma perturbations.
It might be that the field is unstable when $C_{s}^{2}<0$, however, it is
necessary to compute the evolution of the field perturbations and the time of
the stability condition to conclude something in this respect.

On the other hand, the velocity of the field (the time evolution of the field)
has a solution in terms on the scale factor given by:

$v=v_{0}\cdot(1+z)\left(\frac{\left(\frac{1+z_{d}}{1+z}\right)^{m}+1}{\left(1+z_{d}\right)^{m}+1}\right)^{4/m-1}.$
(32)

where $v_{0}$ is the field velocity today.

Finally, the non-canonical term of the action results:

$Q=Q_{0}\cdot(1+z)^{4}\left(\frac{\left(\frac{1+z_{d}}{1+z}\right)^{m}+1}{\left(1+z_{d}\right)^{m}+1}\right)^{4/m+1}\left(\frac{\left(1+z_{d}\right)^{m}-1/3}{\left(\frac{1+z_{d}}{1+z}\right)^{m}-1/3}\right).$
(33)

That is the predicted behaviour for the non-canonical term of the lagrangian
in terms on $z$. More interesting is that $Q(v)$ evolves as radiation in the
first attractor and then acquires a more complicated dependence during the
transition between the attractors, ensuring the continuity of the function.

##  IV Standard cosmological proves

### Age of the universe with this model

The age of the universe according to our model is given by the expression:

$t_{0}=\frac{1}{H_{0}}\int_{0}^{1}\frac{a^{-1}da}{(0.31a^{-3}+0.69f(a))^{1/2}},\\\
t_{0}=1.2987\times 10^{10}\text{years}.$

This is just an approximation of the age of the universe; however, this is an
excellent result taking into account that the parametrization was figured out
for $z<10^{15}$.

On the other hand, the existence of this scalar field implies that the
universe evolves faster than in the standard model, because for a more
negative $\omega$, more accelerated is the expansion and older is the universe
for a given $H_{0}$.

### Matter inhomogeneities evolution

When the Meszaros equation is resolved within the model with the K-essence
scalar field, one solution for the inhomogeneous modes grow with the Hubble
factor de Hubble. Those are the modes which maintain their amplitude after
cross the horizon.

There are many effects on the density of perturbations: it is suppressed the
linear growth $\omega(z)\leq-1/2$, with respect $\Lambda$CDM model, where the
growth is proportional to the scale factor $a$; this suppression raises with
higher $\omega$ at matter domination epoch and entails an earlier beginning to
the dark energy domination epoch (accelerated expansion).

Otherwise, considering the CMB anisotropies spectrum, it is possible to verify
that for values of $m>1.0$, there is formation of the first two acoustics
peaks, then the model has a strong influence during radiation and matter
domination epochs, which not correspond with the observations. On the other
hand, for $m$ values lower than $0.5$, the spectrum goes faster to
$\Lambda$CMB, and they are rejected because do not scale as radiation in that
domination but tend $\omega_{\phi}\sim-1$ for $a<<a_{*}$.

The acoustic peaks correspond to the modes which at the CMB decoupling epoch
were in the maximum compression (odds peaks) or rarefaction (even) and their
position is really sensitive with the state equation of the dark energy.
Actually, the first peak depends on $\Omega_{de}$ monotonically and taking
into account that our model predicts a lower value for this parameter compared
with the $\Lambda$CDM, hence the first peak undergoes a shifting to lower
multipole moments. However, there is a degeneracy that must be broken down in
this phenomenon: the existence of this kind of model implies a rise in the
number of baryons $\Omega_{B}$, that also entails a shifting to lower
multipole moments of the peaks. Furthermore, the larger $\Omega_{B}$ is, the
higher would be the first peak.

Therefore, if $\Omega_{m}$ increases its value within the model, the relative
distance between the peaks decreases, because the mass associated with the
baryons makes that the oscillations occur faster.

###  Statefinder parameters

In order to distinguish between our model and $\Lambda$CDM, rs  has proposed a
test using Statefinder parameters, defined by:

$\displaystyle r$
$\displaystyle=1+\frac{9}{2}\Omega_{\phi}\omega_{\phi}(1+\omega_{\phi})-\frac{3}{2}\Omega_{\phi}\frac{\dot{\omega_{\phi}}}{H},$
(34) $\displaystyle s$
$\displaystyle=1+\omega_{\phi}-\frac{1}{3}\frac{\dot{\omega_{\phi}}}{H\omega_{\phi}},$
(35)

When the parameters are computed, we concluded that the K-essence scalar field
model is slightly off $\Lambda$CDM, because the prediction of the values
$\Omega_{m0}$, $\Omega_{\phi 0}$ and $z_{*}$ are upper than the obtained with
cosmological constant. However, as we have argued in this section, the model
is relaxed and tends to $\Lambda$CDM for redshifts during De-Sitter domination
epoch.

### CMB shift parameter $R$

This parameter measures the shifting of the acoustic peaks from from BAO and
it is defined as the comoving distance between the last scattering surface and
today:

$R=(\Omega_{m}H_{0}^{2})^{1/2}\int_{0}^{1089}\frac{dz}{H(z)},$ (36)

with $H(z)$ depends strongly on the model.
The measured value for this parameter is $R=1.719\pm 0.019$ panotopoulos ,
meanwhile the numerical calculations made with our model give $R_{cal}=1.75$;
the convergence to this value is insured for $z<200$.

##  V Constrains with BBN abundances

As it has been proposed in the introduction of this paper, the main goal with
the effective parametrization of the state equation for the K-essence scalar
field is to avoid the need to fine tuning the initial conditions, but also
reproduce the BBN assess enhancing the predicted abundances. The devise is
really simple, we include the field energy density as relativistic degrees of
freedom during radiation epoch, then compute the Hubble factor and
consequently the Boltzmann equations for the light nuclei.

Firstly, we set out a fiducial baryon to photon ratio $\eta=6.2\times 10^{-9}$
222lambda.gsfc.nasa.gov/product/map/dr4/parameters.cfm and calculate the
capture temperature that deuterium began to be formed via nuclear reactions
given by:

$T_{c,\gamma}=\epsilon_{D}/26=0.088MeV.$ (37)

For this temperature the deuterium is not longer photodissociate bernstein2 ,
where $\epsilon_{D}$ is the deuterium binding energy. Temperatures lower than
$T_{c;\gamma}$ the whole chains related with primordial nucleosynthesis could
be run.

Following the discussion of Bernstein , we assess analytically the value of
the mass fraction $Y_{{}^{4}He}$ and then, using FastBBN kawano ; kellogg ,
Public BBN and BBNreactions bbnreactions , we compute numerically the light
abundances.

In order to introduce the effective degrees of freedom of the K-essence scalar
field, we impose the condition mentioned in section III to restrict the value
of the m-parameter: we want to achieve the largest contribution of the energy
density, been subdominant with respect to the radiation energy density. The
condition can be quantified as:

$\rho_{\phi}|_{rad}=b\cdot\rho_{rad}\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>0\leq b<1.$
(38)

We introduce this energy density in the Hubble parameter during radiation era
and executed the time capture, the neutron’s fraction at this time and
finally, the mass fraction, obtaining the following results:

| $b=0$ | Model with a contribution of $b=0.2$
---|---|---
$H(b)$ ($s^{-1}$) | 1.13 | 1.2379
$t_{c}$ ($s$) | 182 | 169.8
$X_{n}(t_{c})$ | 0.123 | 0.125
$Y_{4}$ | 0.247 | 0.249

The values in the V are compared with the calculated values within the
standard model (where $b=0$, i.e. null contribution of the K-essence scalar
field).

The abundances for the other nuclei are not precisely computed to be reported
using this analytical method, because they depend on the coupled Boltzmann
chains.

However, it is notable that the predicted values for the ${}^{4}He$ is in
agreed with the observational boundaries, therefore our model is an excellent
candidate for dynamical dark energy model in FRWL.

Using FastBBN kawano ; kellogg , Public BBN and BBNreactions bbnreactions , we
compute numerically the light abundances including the relativistic degrees of
freedom of the field as $b$:

Códe | $b$ | $D/H\times 10^{-5}$ | ${}^{3}He/H\times 10^{-5}$ | $Y_{p}$ | ${}^{7}Li/H\times 10^{-10}$
---|---|---|---|---|---
Fast BBN | $0.0$ | $2.335$ | $1.546$ | $0.241$ | $1.268$
Fast BBN | $5.0\times 10^{-3}$ | $2.345$ | $1.548$ | $0.241$ | $1.261$
Fast BBN | $5.0\times 10^{-2}$ | $2.436$ | $1.570$ | $0.245$ | $1.263$
Fast BBN | $0.1$ | $2.537$ | $1.594$ | $0.249$ | $1.261$
Fast BBN | $0.2$ | $2.741$ | $1.639$ | $0.256$ | $1.269$
BBN reactions | $0.0$ | $4.043$ | $2.363$ | $0.243$ | $1.543$
Public Big bang | $0.0$ | $1.542$ | $3.000$ | $0.242$ | $4.884$
Public Big bang | $5.0\times 10^{-3}$ | $1.547$ | $3.002$ | $0.242$ | $4.907$
Public Big bang | $5.0\times 10^{-2}$ | $1.5845$ | $3.015$ | $0.246$ | $5.111$
Public Big bang | $0.1$ | $1.629$ | $3.030$ | $0.250$ | $5.345$
Public Big bang | $0.2$ | $1.717$ | $3.059$ | $0.257$ | $5.825$
our EDE model | $0.2$ | – | – | $0.249$ | –
wmap1  | $0.0$ | $2.75\pm 0.24$ | $0.93\pm 0.055$ | $0.2484\pm 0.0004$ | $3.82\pm 0.66$
wmap2  | $0.0$ | $2.60\pm 0.18$ | $1.04\pm 0.04$ | $0.2479\pm 0.0004$ | $4.15\pm 0.47$

Moreover, it is possible to assess the baryon to photon ratio from the CMB
temperature anisotropies spectrum, taking into account the relation:

$\eta_{B}=\frac{n_{B}}{n_{\gamma}}=5.5\times
10^{-10}\left(\frac{\Omega_{B}h^{2}}{0.022}\right)$ (39)

In the table 4 are shown the values of the parameters related with BBN
($\Omega_{B}$ and $\eta_{B}$) in terms of $b$, where $b=0.0$ is the associated
value to $\Lambda$CDM and $b=0.2$, the maximum contribution of the field
during radiation domination epoch:

$\rho_{\phi}$ | $\Omega_{B}$ | $\Omega_{B}h^{2}$ | $\eta_{B}\times 10^{-10}$
---|---|---|---
$b=0.0$ | 0.044 | 0.02218 | 6.20
$b=0.2$ | 0.053 | 0.02692 | 6.73
wmap1  | $0.0449\pm 0.0028$ | $0.02258\pm 0.00057$ | $6.190\pm 0.145$

According to the results shown in the table V, we plot the abundances for the
light nuclei including ${}^{4}He$, for different contribution of the field,
implemented as effective relativistic degrees of freedom in the code Public
BBN.

Figure 6: $b=0.0$ (solid line), $b=5.0\times 10^{-3}$ (dashed line),
$b=5.0\times 10^{-2}$ (dotted line), $b=0.1$ (width line) and $b=0.2$ (dashed
dotted line).

##  VI Conclusions

In this work, we have made a general description of the dark energy component
as a K-essence scalar field that evolves during radiation as a tracker
($\omega_{\phi}|_{rad}=\frac{1}{3}$) and at late times achieves the De-Sitter
attractor $\omega_{\phi}\rightarrow-1$, emulating cosmological constant
$\Lambda$. All the description was made in the hot cosmological plasm
including matter-radiation components and assuming a spacially flat universe.
With these assumptions for the field, we have proposed an effective
parametrization for the state equation $\omega_{\phi}$ (22), that its free
parameters estimations were obtained by minimizing the function $\chi^{2}$
with the distances modulus of the Type Ia Supernovae.

We have rewritten the velocity of the field and the non-canonical kinetic term
of the K-essence lagrangian in terms of $z$, obtaining the completed behaviour
of the field during the thermal history and this let us say that K-essence
system includes different classes of Quintessence, when it is considered
specific cases of the non-canonical term. Another advantage with this
formulation is that avoids the Fine-tuning of the initials conditions of the
field and its velocity.
As it was expected, during radiation and matter domination epochs the field
has upper predictions for some observables $\Omega_{m0}$, $\Omega_{\phi 0}$
and $z_{*}$ (but respecting the observational bounds for them) comparing with
$\Lambda$CDM model. However, at late times $z\sim 0$ tends asymptotically the
standard model, after it has evolved from the radiation to the second
attractor. All these results are in agreement with the conditions that have
been imposed to resolve the dynamical system kang .
On the other hand, when it is included the non-null contribution of the field
during radiation domination era, the Hubble factor is affected, but also the
time capture and therefore, the primordial light nuclei abundances, because
there were more neutrons out of the equilibrium to form 4He by two body-
reactions. Actually, the whole reactions occurred faster, such that the
production of the nuclei are more effective and drives in an upper mass
fraction.

The predicted value for the 4He abundance prediction according to our model is
inside the observational bounds, then, our model is an excellent candidate to
a dynamic dark energy with a subdominant contribution during radiation and
matter epoch.
Finally, it is remarkable that the results of this paper can be compared with
other kind of models, because the field degrees of freedom can be treated as
effective degrees of some other component (for instance, a Quintessence scalar
field or even, sterile neutrinos lua ). In fact, we have include the field
contribution in the numerical codes in these way, and it let us to compare
degenerations of different kind of models.

###### Acknowledgements.

This work is supported and developed by the Observatorio Astronómico Nacional
under the auspices of Universidad Nacional de Colombia.
We also want to thank Daniel Molano, Carlos Cedeño and all the professors who
give us an advise to enhance this paper in the international events that this
work was presented.

## References

  * (1) S. M. Carroll, Living Rev. Relativity 3 (2001).
  * (2) R. Alpher, H. Bethe, G. Gamov Physics Review 73, 803,(1948).
  * (3) G. Gamov. Physics Review 70, 572,(1946).
  * (4) J. Matsumoto,S. Nojiri. arXiv:hep-th/1001.0220v1,(2010).
  * (5) L. García, J. Tejeiro, L. Castañeda. Proceedings of the International School of Physics Enrico Fermi 178, 309,(2011).
  * (6) A. Rendall. arXiv:gr-qc/0511.158,(2005).
  * (7) C. Armendariz-Picon, V. Mukhanov, P. Steinhardt. arXiv:astro-ph/0006373,(2000).
  * (8) C. Armendariz-Picon, V. Mukhanov, T. Damour. Phys. Lett.B 458, 209),(1999).
  * (9) L. Kawano. Fermilab preprint FERMILAB-PUB-92/04-A,(1992).
  * (10) R. J. Scherrer. Phys. Rev. Lett. 93 011301,(2004).
  * (11) S. Dutta, R. J. Scherrer. arXiv:astro-ph/1006.4166,(2010).
  * (12) R. Putter, E. Linder. arXiv:astro-ph/0705.0400,(2007).
  * (13) R. Saitou, S. Nojiri. arXiv:hep-th/1104.0558,(2011).
  * (14) X. Gao, R. Yang. arXiv:gr-qc/1003.2786v1,(2010).
  * (15) R. Yang, X. Gao. arXiv:gr-qc/1006.4986,(2010).
  * (16) J. Kang, V. Vanchurin, S. Winitzki. Phys. Rev. D 76, 083511,(2007).
  * (17) C. Kellogg. PRadiation Lab preprint OAP-714.
  * (18) J. MacDonald, D. J. Mullan. Physics Rev. D 80, 043507,(2009).
  * (19) A. Riess, A. Filippenko, W. Li, B. Schmidt. Astrophys.J. 536, 62,(2000).
  * (20) R. H. Cyburt, B. D. Fields, K. A. Olive. Physics Letters B 567, 227-234,(2003).
  * (21) A. Coc, E. Vangioni-Flam, P. Descouvemont, A. Adahchour, C. Angulo. The Astrophysical Journal, 600:544 552,(2001).
  * (22) E. Komatsu, et. al. ApJS, 192, 18,(2010).
  * (23) E. J. Copeland, M. Sami, S. Tsujikawa. arXiv:hep-th/0603057 ,(2006).
  * (24) G. Steigman. Proceedings IAU Symposium No. 265, (2009).
  * (25) J. Bernstein, L. S. Brown, G. Feinberg. Reviews of Modern physics, Vol. 61, No. 1,(1989).
  * (26) G. Steigman, J.P. Kneller. arXiv:astro-ph/0210500,(2002).
  * (27) G. Panotopoulos. arXiv:gr-qc/1107.4475v1,(2011).
  * (28) T. Chiba. arXiv:astro-ph/0206298,(2002).


# Graphical Representations for Algebraic Constraints
of Linear Structural Equations Models

Thijs van Ommen<EMAIL_ADDRESS>
Information and Computing Sciences, Utrecht University, The Netherlands
Mathias Drton<EMAIL_ADDRESS>
Mathematics, Technical University of Munich, Germany

###### Abstract

The observational characteristics of a linear structural equation model can be
effectively described by polynomial constraints on the observed covariance
matrix. However, these polynomials can be exponentially large, making them
impractical for many purposes. In this paper, we present a graphical notation
for many of these polynomial constraints. The expressive power of this
notation is investigated both theoretically and empirically.

Keywords: structural equation model; mixed graph; multivariate Gaussian;
algebraic constraint.

##  1 Introduction

It is well known that for any directed acyclic graph, the set of observational
distributions realizable by that graph can be described in terms of a set of
conditional independence constraints that are imposed on the distributions.
Our understanding of these constraints has enabled many sophisticated
algorithms for learning the structure of Bayesian networks.

When latent confounders are considered, as is important for purposes of causal
discovery, the situation changes and additional constraints such as the Verma
constraint (Robins, 1986; Verma and Pearl, 1991) are needed in addition to
conditional independences. Here we need to distinguish between the
nonparametric (as well as discrete) case on the one hand and the linear case
on the other. In the nonparametric case, the constraints involved are becoming
better and better understood (Tian and Pearl, 2002; Shpitser et al., 2014). In
the linear case, specifically for linear structural equation models (LSEMs),
these same constraints apply (Shpitser et al., 2018), as well as others such
as the vanishing tetrad constraint (Spirtes et al., 2000) and its higher-
dimensional generalizations (Sullivant et al., 2010).

In principle, the constraints imposed by a graph in the linear case can be
found using methods from computer algebra (Drton, 2018), yet these can be
extremely slow in practice. Further, the constraints are output in the form of
long polynomials, which do not immediately lend themselves to a better
understanding of these graphical models. In this paper, we present a new,
graphical description of such constraints, which is much more succinct and
which we believe will help further the understanding of linear structural
equation models. In particular, knowledge about constraints helps checking
whether different graphs encode observationally indistinguishable hypotheses
because the graphs are equivalent in terms of defining the same statistical
model. Moreover, constraints may be deployed for goodness-of-fit tests (e.g.,
Leung and Drton (2018)), and they can also lead to new criteria to decide
parameter identifiability (Weihs et al., 2017).

This paper is structured as follows. Section 2 discusses preliminaries about
linear structural equation models. In Section 3, we state our contributions:
the definition of graphical constraints, an algorithm for deriving a list of
constraints for a given model, and a number of ancillary results. Empirical
results are given in Section 4, and Section 5 concludes.

##  2 Preliminaries

In a linear structural equation model (LSEM), the value of each variable
$X_{v}$ is governed by an equation of the form (see e.g. Foygel et al. (2012)
for a more extensive exposition)

$X_{v}=\lambda_{0,v}+\sum_{w}\lambda_{w,v}X_{w}+\epsilon_{v}.$

The structure of an LSEM can be represented by a _mixed graph_ $G=(V,D,B)$,
where the nodes $V$ index the random variables $X_{v}$ of the LSEM, each
directed edge $v\to w\in D$ denotes that $X_{v}$ appears on the right-hand
side of $X_{w}$’s structural equation, and $v\leftrightarrow w\in B$ denotes
that the noise terms $\epsilon_{v}$ and $\epsilon_{w}$ may have nonzero
correlation. Note that LSEMs may contain directed cycles, though self-loops
($v\to v$) are not allowed. Also allowed are _bows_ : the co-occurrence of a
directed and a bidirected edge between the same pair of nodes.

Ignoring the (to us) irrelevant additive terms $\lambda_{0,v}$, the parameters
of an LSEM are the matrices $\Lambda$ and $\Omega$. $\Lambda$ must have
$\Lambda_{v,w}=0$ wherever $v\to w\notin D$, and must be such that $I-\Lambda$
is invertible. This last condition is automatically satisfied if $G$ contains
no directed cycles. We write $\mathbb{R}^{D}_{\text{reg}}$ for the space of
such $\Lambda$’s. $\Omega$ must have $\Omega_{v,w}=0$ except on the diagonal
and when $v\leftrightarrow w\in B$, and further it must be positive definite.
We write $\mathrm{PD}(B)$ for the space of such $\Omega$’s, and write
$\mathrm{PD}_{n}$ for the space of $n\times n$ positive definite matrices.

After observing samples from an LSEM, we can estimate the covariance matrix
$\Sigma$ on the observed variables $(X_{v})_{v\in V}$. Two central problems
ask what we can do when we know $\Sigma$: the _(parameter) identification_
problem is to reconstruct the parameters $\Lambda$ and $\Omega$ from $\Sigma$
when the graph is known, and the _structure learning_ problem is to find an
appropriate graph when only $\Sigma$ is known.

For structure learning, it is important to know for each graph $G$ what set of
$\Sigma$’s could be observed for that graph and any choice of parameters. This
set, the _model_ of $G$, is denoted $\mathcal{M}(G)$. In general,
$\mathcal{M}(G)$ is a semialgebraic set (we refer to Cox et al. (2015) for
further background on algebraic geometry), described by polynomial equalities
and inequalities. We will only study the equalities ($f(\Sigma)=0$) in this
paper; as argued by van Ommen and Mooij (2017), the inequalities are of
secondary importance for structure learning. The set of polynomials that are
zero everywhere on $\mathcal{M}(G)$ form an _ideal_ $I_{G}$. This ideal can be
described by a finite number of _generators_ (our constraints play this role).
The _algebraic model_ $\overline{\mathcal{M}(G)}$ is defined as $V(I_{G})$,
where $V(I)$ denotes the set of points on which all polynomials in ideal $I$
evaluate to zero. In other words,
$\overline{\mathcal{M}(G)}\supseteq\mathcal{M}(G)$ ignores any inequality
constraints imposed by $\mathcal{M}(G)$.

Each of the constraints describing an algebraic model is an _homogeneous_
polynomial: each term has the same total degree. In fact, these polynomials
obey a stronger property we call _$V$ -homogeneity_: for each model variable
$v\in V$, the number of occurrences of $v$ in the indices of the $\sigma$’s
(counting $\sigma_{vv}$ twice) is the same for all terms in the polynomial.
Rescaling an LSEM variable by a constant does not affect whether $f(\Sigma)=0$
iff $f$ is $V$-homogeneous.

###  2.1 Rational and Polynomial Constraints

In order to derive polynomial constraints, we will need expressions for the
parameters as a function of $\Sigma$. This is the problem of parameter
identification for LSEMs, and a powerful method for this problem is the half-
trek criterion (HTC: Foygel et al., 2012). A _half-trek_ is a directed path,
except that the first edge on the path may instead be bidirected. The set of
nodes reachable from $v$ by a half-trek is denoted $\mathrm{htr}(v)$. The HTC-
identifiability theorem states that if for each $v\in V$, a set $Y_{v}$ exists
that satisfies certain conditions, then rational functions exists that
identify the parameters. We write $\mathcal{Y}$ for the sequence $(Y_{v})_{v}$
of HTC-identifying sets, and $\Lambda_{\mathcal{Y}}(\Sigma)$ for the function
mapping $\Sigma$ to $\Lambda$.

We build on Theorem 1 of van Ommen and Mooij (2017), which shows that for an
HTC-identifiable graph $G$ with HTC-identifying sets
$\mathcal{Y}=(Y_{v})_{v}$:

  1. 1.

For generic $\Sigma\in\mathcal{M}(G)$, $\Lambda_{\mathcal{Y}}(\Sigma)$ is
defined, is in $\mathbb{R}^{D}_{\text{reg}}$, and satisfies the _rational
constraints_

$[(I-\Lambda_{\mathcal{Y}}(\Sigma))^{T}\Sigma(I-\Lambda_{\mathcal{Y}}(\Sigma))]_{v,w}=0$
(1)

for all $v\neq w$ with $v\leftrightarrow w\notin B$, $v\notin Y_{w}$, and
$w\notin Y_{v}$;

  2. 2.

All $\Sigma\in\mathcal{M}(G)$ satisfy the _polynomial constraints_ obtained by
multiplying out the denominators from (1);

  3. 3.

Among $\Sigma$ that satisfy the rational constraints (1) and have
$\Lambda_{\mathcal{Y}}(\Sigma)\in\mathbb{R}^{D}_{\text{reg}}$, generically
$\Sigma\in\mathcal{M}(G)$.

Statements 1 and 2 (showing that $\Sigma\in\mathcal{M}(G)$ implies
satisfaction of constraints) are strong: statement 2, about the polynomial
constraints, is not even limited to generic $\Sigma$. Statement 3 in the other
direction is weaker. In particular, it does not imply anything for $\Sigma$
where the rational constraints are undefined, or for which
$\Lambda_{\mathcal{Y}}(\Sigma)\notin\mathbb{R}^{D}_{\text{reg}}$. We will
revisit this issue in Section 3.3.

##  3 Main Results

In this section, we present the definition of graphical constraints, followed
by an algorithm that finds such constraints for any HTC-identifiable mixed
graph. Section 3.3 discusses the question of how well a list of graphical
constraints describes an algebraic model, and Section 3.4 introduces
transformation operations that may be performed on graphical constraints.

###  3.1 Graphical Constraints

A _graphical constraint_ is an undirected bipartite graph (usually a tree), in
which each node is labelled with a subset of the model variables $V$.

Each graphical constraint represents a polynomial. This polynomial is the
symbolic determinant of the matrix $M$ defined as follows. Let $A$, $B$ be the
two parts of the bipartite graphical constraint. The matrix $M$ has a row for
each pair $(a,v)$ with $a\in A$ and $v$ in $a$’s label, and a column for each
pair $(b,w)$ with $b\in B$ and $w$ in $b$’s label. We use these pairs to index
the rows and columns of $M$. We require that $M$ is square. Then $M$ is given
by

$M_{(a,v),(b,w)}=\begin{cases}\sigma_{v,w}&\text{if $a$ and $b$ are
adjacent;}\\\ 0&\text{otherwise.}\end{cases}$ (2)

We say an observational covariance matrix $\Sigma$ _satisfies_ a graphical
constraint if the determinant of $M$ evaluates to 0 at $\Sigma$.

A different ordering of the rows and columns of $M$ may flip the sign of the
determinant. Because a graphical constraint does not fix a particular
ordering, it only represents a polynomial up to sign. This does not pose a
problem, since it has no effect on whether or not for given $\Sigma$, the
polynomial equals 0.

acbc

(a) $\begin{bmatrix}\sigma_{ab}&\sigma_{ac}\\\
\sigma_{cb}&\sigma_{cc}\end{bmatrix}$ (b) $\mathstrut a$$\mathstrut
b$$\mathstrut c$$\mathstrut d$(c) dabacd(d)
$\begin{bmatrix}\sigma_{da}&\sigma_{dc}&0\\\
\sigma_{aa}&\sigma_{ac}&\sigma_{ad}\\\
\sigma_{ba}&\sigma_{bc}&\sigma_{bd}\end{bmatrix}$ (e)

Figure 1: (a) a graphical constraint (we draw these in a zig-zag way to make
the bipartite parts more easily recognizable); (b) the matrix represented by
this graphical constraint; (c) an LSEM (with a bow and a directed cycle)
imposing a single constraint, namely the graphical constraint displayed in
(d), with (e) the matrix it represents.

Two examples of graphical constraints are displayed in Figure 1 (a) and (d).
The determinant of (b) is $\sigma_{ab}\sigma_{cc}-\sigma_{ac}\sigma_{bc}$,
which is zero iff the partial correlation of $a$ and $b$ given $c$ is zero;
this is the Gaussian analogue of conditional independence.

We say a graphical constraint is in normal form if no two nodes have the same
set of neighbours. A pair of nodes that have the same set of neighbours can be
merged into a single node whose label is the union of the original nodes’
labels, without changing the matrix. Restricting graphical constraints to
normal form ensures that two different graphical constraints do not describe
the same matrix. For example, when the entire tree consists of a central node
with some leaf nodes around it, the matrix (2) will have no zero entries, and
many similar trees represent the same matrix.

The polynomial represented by a graphical constraint may have exponentially
many terms. An immediate advantage of the graphical representation is that
evaluating a polynomial for given $\Sigma$ can be made computationally much
more efficient by using algorithms for computing determinants (of sparse
matrices), compared to evaluating the exponentially many terms. This is useful
for testing a constraint on observational data.

###  3.2 Construction Algorithm

Algorithm 1 can be used to find a graphical representation for each of the
constraints in (1). The relation between rational/polynomial constraints and
the graphical constraints is made precise below in Theorem 1.

Input: An HTC-identifiable graph $G$, a family $\mathcal{Y}=(Y_{v})_{v}$ of
HTC-identifying sets, and a pair $\\{v,w\\}$ as in (1)

Output: A graphical constraint

1 Start with a graph consisting of two ‘seed’ nodes $t_{1}$ and $t_{2}$ joined
by an edge, one with label $\\{v\\}$ and one with label $\\{w\\}$;

2 expand(_$t_{1}$_);

3 expand(_$t_{2}$_);

4 Wherever two or more leaf nodes are attached to the same node, merge them
into a single leaf node labelled with the union of those leaf nodes’ labels.
Subroutine  _expand(_$t$_)_:

5 Let $\\{v\\}$ be the label of $t$;

6 Replace the label of $t$ by $\\{v\\}\cup\mathrm{pa}(v)$;

7 For each $y\in Y_{v}$, create a new node $t_{y}$ attached to $t$ with label
$\\{y\\}$;

8 For each $y\in Y_{v}$ with $y\in\mathrm{htr}(v)$, expand(_$t_{y}$_);

end

Algorithm 1 HTC-based graphical constraint construction.

Line 4 of the algorithm brings the constraint into normal form. It suffices to
check leaf nodes because the graph is a tree, so the only pairs of nodes that
could have the same set of neighbours are leaf nodes.

Figure 2 displays an example where Algorithm 1 is applied to the graph (a).
This graph is bow-free and acyclic, in which case it is possible to take
$Y_{v}=\mathrm{pa}(v)$ as the HTC-identifying sets $\mathcal{Y}$. With this
input, the algorithm produces the graphical constraint in (b). This constraint
represents the matrix (c), whose determinant is a polynomial of degree 3 with
4 terms. The LSEM in Figure 1 (c) is not bow-free or acyclic; with
$Y_{a}=\\{c,d\\}$, $Y_{b}=Y_{c}=\\{d\\}$ and $Y_{d}=\emptyset$, the algorithm
outputs the graphical constraint in Figure 1 (d).

$\mathstrut a$$\mathstrut b$$\mathstrut c$$\mathstrut d$

(a) ⓒbⓓaba(b) $\begin{bmatrix}\sigma_{cb}&\sigma_{cd}&0\\\
\sigma_{ab}&\sigma_{ad}&\sigma_{aa}\\\
\sigma_{bb}&\sigma_{bd}&\sigma_{ba}\end{bmatrix}$ (c)

Figure 2: (a) an LSEM; (b) the graphical constraint output by Algorithm 1
(with the original seed node labels marked by circles); (c) the matrix
represented by this graphical constraint.

###### Theorem 1

The polynomial represented by the graphical constraint from Algorithm 1 is
either equal to or a multiple of the numerator of the rational constraint (1).

Proof  At each point in the algorithm _after_ an expansion step, we have a
tree graph in which some nodes might still have to be expanded. Take the seed
node that is expanded first to be the root of each of these trees. To each of
these intermediate trees, we associate a matrix like (2), except that at non-
zero positions $y,v$ or $v,y$ where $y$ is to-be-expanded, the matrix entry is
not $\Sigma_{y,v}$ but $[(I-\Lambda_{\mathcal{Y}}(\Sigma))^{T}\Sigma]_{y,v}$.
Here and in the following, we make slight abuse of notation and take $v$,
$\mathrm{pa}(v)$, $Y_{v}$, etc. to refer both to labels (i.e. elements of $V$)
and to the specific row or column indices corresponding to the (node, label
element)-pair under consideration.

The proof proceeds by induction. Our induction hypothesis is that after each
expansion step in the algorithm, the matrix at that point has a determinant
that is a multiple of (1).

First we prove that for any vector
$\mathbf{c}=[c_{v},c_{p_{1}},\ldots,c_{p_{n}}]^{T}$ where
$\\{p_{1},\ldots,p_{n}\\}=\mathrm{pa}(v)$,

$\lvert\mathbf{A}^{(v)}\rvert\cdot[(I-\Lambda_{\mathcal{Y}}(\Sigma))_{\\{v\\}\cup\mathrm{pa}(v);v}]^{T}\mathbf{c}=\begin{vmatrix}c_{v}&c_{p_{1}}\cdots
c_{p_{n}}\\\ \mathbf{b}^{(v)}&\mathbf{A}^{(v)}\end{vmatrix}.$ (3)

Here $\mathbf{A}^{(v)}$ and $\mathbf{b}^{(v)}$ are the matrix and vector from
(Foygel et al., 2012, proof of Theorem 1), where they are used to identify
$\Lambda_{\mathcal{Y}}(\Sigma)_{\mathrm{pa}(v),v}$ as the solution of
$\mathbf{A}^{(v)}\cdot\Lambda_{\mathcal{Y}}(\Sigma)_{\mathrm{pa}(v),v}=\mathbf{b}^{(v)}$.
They are defined by $\mathbf{A}^{(v)}=\Sigma_{Y_{v},\mathrm{pa}(v)}$ and
$\mathbf{b}^{(v)}=\Sigma_{Y_{v},v}$, except in rows where
$y\in\mathrm{htr}(v)$: there
$\mathbf{A}^{(v)}_{y,\mathrm{pa}(v)}=[(I-\Lambda_{\mathcal{Y}}(\Sigma))^{T}\Sigma]_{y,\mathrm{pa}(v)}$
and
$\mathbf{b}^{(v)}_{y}=[(I-\Lambda_{\mathcal{Y}}(\Sigma))^{T}\Sigma]_{y,v}$.
For each $p\in\mathrm{pa}(v)$, by Cramer’s rule,
$\Lambda_{\mathcal{Y}}(\Sigma)_{p,v}=\lvert\mathbf{A}^{(v)}_{p}\rvert/\lvert\mathbf{A}^{(v)}\rvert$,
where $\mathbf{A}^{(v)}_{p}$ is the matrix obtained from $\mathbf{A}^{(v)}$ by
replacing column $p$ by $\mathbf{b}^{(v)}$. So

$\displaystyle\lvert\mathbf{A}^{(v)}\rvert\cdot[I-\Lambda_{\mathcal{Y}}(\Sigma)]_{\\{v\\}\cup\mathrm{pa}(v);v}$
$\displaystyle=\left[\lvert\mathbf{A}^{(v)}\rvert,-\lvert\mathbf{A}^{(v)}_{p_{1}}\rvert,\ldots,-\lvert\mathbf{A}^{(v)}_{p_{n}}\rvert\right]^{T}.$

For $\mathbf{c}$ a standard basis vector, (3) can be verified by noting that
for $\mathbf{c}$ with $c_{p}=1$ ($p\in\mathrm{pa}(v)$),

$-\lvert\mathbf{A}^{(v)}_{p}\rvert=-\begin{vmatrix}1&\mathbf{0}\\\
\mathbf{0}&\mathbf{A}^{(v)}_{p}\end{vmatrix}=\begin{vmatrix}0&c_{p_{1}}\cdots
c_{p_{n}}\\\ \mathbf{b}^{(v)}&\mathbf{A}^{(v)}\end{vmatrix},$

where the sign flips due to the transposition of two columns. Now the case of
general $\mathbf{c}$ follows by multilinearity of the determinant.

_(base case)_ The first expansion performed by the algorithm turns a tree with
only the two seed nodes $v$ and $w$ into a tree where only root node $v$ has
been expanded, so all other nodes in the tree are children of $v$. The
corresponding matrix has columns $\\{v\\}\cup\mathrm{pa}(v)$ and rows
$\\{w\\}\cup Y_{v}$ (note that the roles of rows and columns can be
interchanged). Because $w$ is to-be-expanded, entry $w,x$ in $w$’s row equals
$[\Sigma(I-\Lambda_{\mathcal{Y}}(\Sigma))]_{x,w}$. For each $y\in Y_{v}$, that
$y$ is to-be-expanded iff $y\in\mathrm{htr}(v)$, so the other rows of the
matrix equal either entries of $\Sigma$ or of
$(I-\Lambda_{\mathcal{Y}}(\Sigma))^{T}\Sigma$, the choice agreeing with the
choice in the definitions of $\mathbf{A}^{(v)}$ and $\mathbf{b}^{(v)}$,
establishing that this matrix equals the one in the right-hand side of (3) for
$c_{x}=[\Sigma(I-\Lambda_{\mathcal{Y}}(\Sigma))]_{x,w}$. By (3), the
determinant equals a multiple of

$[(I-\Lambda_{\mathcal{Y}}(\Sigma))^{T}\Sigma(I-\Lambda_{\mathcal{Y}}(\Sigma))]_{v,w},$

proving that the induction hypothesis holds after the first expansion.

_(induction step)_ For each subsequent expansion, let $v$ be the node that is
expanded in this step, so that its label after expansion is
$\\{v\\}\cup\mathrm{pa}(v)$ and the union of its children’s labels is $Y_{v}$.
Write $L$ for the label of $v$’s parent in the tree (it is the same before and
after the expansion). The matrix $M^{\prime}$ before expansion differs from
the matrix $M$ after expansion in that the rows and columns corresponding to
$\mathrm{pa}(v)$ and $Y_{v}$ are absent in $M^{\prime}$, and the entries
$v,\ell$ (or $\ell,v$) for $\ell\in L$ are equal to
$[(I-\Lambda_{\mathcal{Y}}(\Sigma))^{T}\Sigma]_{v,\ell}$. We need to show that
the determinant of the matrix $M$ is a multiple of the determinant of the
matrix $M^{\prime}$.

Because $\lvert\\{v\\}\cup\mathrm{pa}(v)\rvert=\lvert Y_{v}\rvert+1$, each
term in the expansion of the determinant of $M$ contains exactly one entry
from $(\\{v\\}\cup\mathrm{pa}(v))\times L$. So

$\displaystyle\lvert M\rvert$ $\displaystyle=\sum_{\ell\in L}\lvert
M_{\text{with $(\\{v\\}\cup\mathrm{pa}(v))\times(L\setminus\\{\ell\\})$
replaced by zeros}}\rvert$ $\displaystyle=\sum_{\ell\in L}\lvert
M_{(\\{v\\}\cup\mathrm{pa}(v));(\\{\ell\\}\cup Y_{v})}\rvert\cdot\lvert
M_{\text{other rows and columns}}\rvert$
$\displaystyle=\lvert\mathbf{A}^{(v)}\rvert\sum_{\ell\in
L}[(I-\Lambda_{\mathcal{Y}}(\Sigma))_{\\{v\\}\cup\mathrm{pa}(v);v}]^{T}\Sigma_{\\{v\\}\cup\mathrm{pa}(v);\ell}\cdot\lvert
M_{\text{other rows and columns}}\rvert$
$\displaystyle=\lvert\mathbf{A}^{(v)}\rvert\cdot\lvert M^{\prime}\rvert.$

This completes the proof.

This theorem shows that whenever the polynomial constraint from Section 2.1,
point 2 is satisfied, then so is the corresponding graphical constraint.

The proof also makes explicit by what factors the rational constraint is
multiplied: for each expansion of a node with label $\\{v\\}$, there is a
multiplication by $\lvert\mathbf{A}^{(v)}\rvert$. As the final product is a
polynomial, clearly the product of these factors is a multiple of the rational
constraint’s denominator. So if the rational constraints are undefined for
some $\Sigma$, this must mean that a relevant
$\lvert\mathbf{A}^{(v)}\rvert=0$. By inspecting the construction, we can also
isolate from the graphical constraint the part that represents the polynomial
$\lvert\mathbf{A}^{(v)}\rvert$: it is the graphical constraint obtained by
starting with a single node with label $\\{v\\}$ and expanding it, then
finally removing the label $v$ from the starting node.

###  3.3 Describing Algebraic Models Using Graphical Constraints

In terms of algebraic geometry (for which we refer to Cox et al. (2015)), the
algebraic model $\overline{\mathcal{M}(G)}$ is an irreducible variety, which
is described by a prime ideal. The constraints found by Algorithm 1 generate
an ideal that is not necessarily prime. There are two ways in which an ideal
can fail to be prime, using that an ideal is prime iff it is _radical_ and
_primary_ :

  * •

An ideal $I$ is nonradical if $f^{k}\in I$ for some polynomial $f$ and some
$k\geq 2$, but $f\notin I$. The radical $\sqrt{I}$ of the ideal $I$ also
includes all such polynomials $f$. We are mostly concerned with the ideal’s
ability to describe a model, and $V(I)=V(\sqrt{I})$, so failure to be radical
does not impact our goal.

  * •

If an ideal $I$ is nonprimary, then $V(I)$ is the union of multiple
irreducible varieties. For instance, the ideal generated by $xy$ describes a
variety that is the union of the two hyperplanes $x=0$ and $y=0$. Primary
ideals describing these irreducible varieties can be found computationally
using _primary decomposition_ algorithms, but these algorithms can be
extremely computationally expensive. It follows from Theorem 1 of van Ommen
and Mooij (2017) that one of the component varieties is the algebraic model we
want to describe, meaning that the remaining components are spurious.

The spurious components of a nonprimary ideal lead to false positives: sets of
$\Sigma$ that satisfy the polynomial constraints but do not belong to the
algebraic model. We may further classify ideals by the types of spurious
primary components to which they give rise:

  * •

If an ideal $I_{\text{comp}}$ in the primary decomposition contains (a power
of) a principal minor of $\Sigma$, then $V(I_{\text{comp}})$ only consists of
$\Sigma$ for which this principal minor evaluates to 0. Such $\Sigma$ are on
the boundary of the positive definite cone, so
$V(I_{\text{comp}})\cap\mathrm{PD}_{n}=\emptyset$ and this spurious component
does not impact the constraints’ ability to distinguish which covariance
matrices are in the algebraic model. If all spurious components of an ideal
have this property, we call the ideal _PD-primary_.

  * •

Suppose some $f\in I_{\text{comp}}$ contains a monomial
$\prod_{v}\sigma_{vv}^{k_{v}}$ for some sequence $(k_{v})_{v}$; i.e., a
monomial consisting entirely of elements from the diagonal of $\Sigma$. If $f$
is $V$-homogeneous, there can be at most one monomial of this form. Then the
identity matrix $I$, or in fact any diagonal $\Sigma\in\mathrm{PD}_{n}$, is
not in $V(I_{\text{comp}})$, because
$f(\Sigma)=\prod_{v}\sigma_{vv}^{k_{v}}>0$. If all spurious components of an
ideal contain such a monomial, we call the ideal _$I$ -primary_.

All algebraic models we consider in this paper _do_ contain all such diagonal
$\Sigma$ — even the smallest model, described by the empty mixed graph,
contains all these. This makes it clear that this component is spurious. So in
this case, it may be that
$V(I_{\text{comp}})\cap\mathrm{PD}_{n}\neq\emptyset$, but we do have that the
polynomial constraints contain all the information needed to uniquely indicate
an algebraic model.

  * •

In more general cases, the spurious component may be equal to an algebraic
model other than the one we mean to describe. We may be able to distinguish
which component is correct and which are spurious by referencing the
denominators in the rational constraints (1). But the polynomial constraints
by themselves may not contain enough information to indicate a single
algebraic model.

Note that primary ideals are PD-primary, and PD-primary ideals are
$I$-primary.

As mentioned in Section 2.1, Theorem 1 of van Ommen and Mooij (2017) leaves
open the possibility of $\Sigma$ that satisfy the polynomial constraints but
are not in $\mathcal{M}(G)$, when either $\Lambda_{\mathcal{Y}}$ is undefined
at $\Sigma$ or
$\Lambda_{\mathcal{Y}}(\Sigma)\notin\mathbb{R}^{D}_{\text{reg}}$. This final
case can only occur if $G$ contains a directed cycle, and we will not address
it further in this paper. The case of $\Lambda_{\mathcal{Y}}$ being undefined
arises when a matrix $\mathbf{A}^{(v)}$ in the HTC identification procedure is
singular (see the end of Section 3.2). We can use this to show that certain
ideals are well-behaved:

###### Proposition 2

If $G$ is ancestral, the ideal generated by the graphical constraints of
Algorithm 1 with the choice $Y_{v}=\mathrm{pa}(v)$ for each $v$ is PD-primary.

Proof  In an ancestral graph, $\mathrm{pa}(v)\cap\mathrm{htr}(v)=\emptyset$,
so the expand subroutine will not call itself recursively. Then the graphical
constraint representing $\lvert\mathbf{A}^{(v)}\rvert$ (see the end of Section
3.2) has two nodes, both with label $\mathrm{pa}(v)$. This describes a
principal minor, implying that the ideal in question is PD-primary.

###### Proposition 3

If $G$ is bow-free and acyclic, the ideal generated by the graphical
constraints of Algorithm 1 with the choice $Y_{v}=\mathrm{pa}(v)$ for each $v$
is $I$-primary.

Proof  Since $G$ is acyclic, any spurious components must be contained in
regions where the rational constraints are not defined, meaning they must each
lie within one of the varieties defined by $\lvert\mathbf{A}^{(v)}\rvert=0$
for some $v$. For the choice $Y_{v}=\mathrm{pa}(v)$, we see from the graphical
constraint representing $\lvert\mathbf{A}^{(v)}\rvert$ that the expansion of
$\lvert\mathbf{A}^{(v)}\rvert$ contains a unique monomial
$\prod_{v}\sigma_{vv}^{k_{v}}$. This implies that the ideal in question is
$I$-primary.
We will see in Section 4 that these sufficient conditions are far from
necessary.

###  3.4 Transformation of Graphical Constraints

Performing elementary row and column operations on the matrix (2) does not
change its determinant except by a constant factor. In particular, the
operation of adding or subtracting one row or column from another will yield a
matrix that can again be represented by a graphical constraint, if the two
involved rows/columns were labelled with the same LSEM variable and the
entries of the resulting matrix are still all either $0$ or $\sigma_{vw}$.
Then these graphical constraints are equivalent, even though they are not
isomorphic as graphs.

An application of these matrix transformations is to remove spurious factors.
Theorem 1 shows that the polynomial represented by a graphical constraint can
be equal to the numerator of the corresponding rational constraint, or it can
contain additional factors. Getting rid of such factors can sometimes remove
problematic primary components from the ideal. If a matrix can be transformed
into a $2\times 2$ block diagonal matrix, this shows that its determinant
factorizes into the determinants of the two diagonal blocks. Below we
illustrate a class of such steps that can be recognized graphically.

$R$$XA$$XB$$\tau$$\tau$$\lvert A\rvert$

(a) is transformed into $R$$X$$XAB$$\tau$$\tau$$+\lvert A\rvert$(b)

Figure 3: (a) Transformation preconditions: all trees $\tau$ at the right node
have copies at the left node and the weight of the left edge equals $\lvert
A\rvert$; (b) result of the transformation.

This transformation applies to tree-shaped constraints and revolves around
three nodes (called left, right, and central here). Let $X$ be the
intersection of left’s and right’s labels, and $A$ and $B$ the labels in left
but not in right and vice versa. The transformation’s preconditions are
illustrated in Figure 3 (a): For each subtree at the right node (not
containing the central node), there is an identical subtree at the left node;
and the weight of the edge between the central and left nodes equals $\lvert
A\rvert$. Here the _weight_ of an edge counts the number of occurrences in a
term of the determinant of the elements from the submatrix of $M$ with rows
and columns corresponding to the endpoints of this edge; for a tree-shaped
graphical constraint, this number is the same for all terms in the
determinant.

The result of the transformation is displayed in Figure 3 (b): the edge
between the central and left nodes is deleted, and the nonshared labels of the
left node are removed there and added to the right node.

###### Proposition 4

The transformation described above preserves the represented polynomial up to
sign.

Proof  Assume w.l.o.g. that the left and right nodes are columns. Add ‘cross-
edges’ between the left node and the roots of the right $\tau$-subforest by
adding to each row in the right $\tau$-subforest the corresponding row from
the left $\tau$-subforest, and similarly reducing columns in the left
$\tau$-subforest by their right counterparts. Then reduce the columns of the
left node’s shared labels $X$ by their right counterpart. Putting the left
node with all attached parts first in the matrix and the central and right
nodes with attached parts second, we find that the diagonal blocks are square
and that all entries in the bottom left block are 0, so making the top right
block 0 as well does not change the determinant.

##  4 Computational Study

Proposition 3 shows that for bow-free acyclic graphs, with the choice
$Y_{v}=\mathrm{pa}(v)$, Algorithm 1 is guaranteed to produce an $I$-primary
ideal. Since PD-primary ideals are preferable, we ran a computational study
over all bow-free acyclic graphs on five variables with nine edges. (For these
graphs, the model is described by exactly one constraint, and the type of
ideal generated by this constraint can be determined by its factorization.
Graphs with fewer edges impose more constraints, the analysis of which was not
computationally feasible.)

The results were as follows: Without transformation, among the 86 algebraic
equivalence classes (counting isomorphic classes only once, so e.g. the class
imposing graphical constraint $a$—$b$ and the one imposing $d$—$e$ are
together counted once), there were 5 classes containing mixed graphs for which
the algorithm produced a non-PD-primary ideal (each of these classes also
contained graphs producing PD-primary ideals). Using the transformation
displayed in Figure 3, all non-PD-primary ideals were converted to PD-primary
ones. So in practice, the ideals are often of a higher quality than could be
guaranteed by Proposition 3.

The above is a measure of the utility of Algorithm 1 and of the transformation
from Figure 3. We also want to know whether the class of graphical constraints
has enough expressive power to represent arbitrary polynomial constraints,
preferably without spurious factors. To study this, for each algebraic
equivalence class on five nodes with at least nine edges per mixed graph,
allowing bows but not directed cycles, we tried the algorithm on all mixed
graphs in the class and all HTC-identifying choices of $\mathcal{Y}$. If we
did not find a graphical constraint without spurious factors this way, we also
tried various transformations or even brute-force search. We did the same for
four-node graphs with at least five edges allowing bows and directed cycles.

primary | PD-primary | 5-node acyclic | 4-node general
---|---|---|---
tree | tree | 374 | 16
nontree(?) | tree | 2 | 0
none | tree | 2 | 0
? | tree | 3 | 0
? | ? | 17 | 3
total | 398 | 19
Table 1: Number (up to isomorphism) of algebraic equivalence classes imposing
one constraint for which certain types of graphical constraints are known to
exist.

The results are summarized in Table 1. Among the five-node graphs, there were
18 equivalence classes without any HTC-identifiable members, so Algorithm 1
could not be used. One of these classes imposed a polynomial constraint of
degree 7. For this constraint we found a graphical form via brute-force
search. The other 17 equivalence classes impose constraints of degrees 9
through 14, for which such a brute-force search was infeasible.

It is noteworthy that in all classes where we know at least one graphical
constraint, we know a (tree-shaped) PD-primary constraint. In the majority of
classes, this constraint is in fact primary. In two cases, all primary
graphical constraints are nontrees (confirmed by brute-force search in one of
these two cases), and in two other cases, a brute-force search found no
primary graphical constraints of any kind. In three further cases, primary
graphical constraints have not been found but might exist.

##  5 Conclusion

We have presented a graphical formalism to give succinct representations of
polynomials that appear as constraints on the algebraic model of an LSEM, and
illustrated the expressive power of these graphical constraints for purposes
of describing algebraic models.

We also presented an algorithm that derives a list of graphical constraints
for any HTC-identifiable graph, and a transformation that can be used to
simplify some graphical constraints. While our empirical results show they are
quite powerful at their tasks, there is room for improvement here. A number of
generalizations of HTC have been proposed, such as EID and TSID (Weihs et al.,
2017), and generalizing our algorithm to work with these would enable the
derivation of graphical constraints for a larger class of mixed graphs.
Improved algorithms might also avoid some of the spurious factors which are
sometimes produced by our current algorithm.

Another important direction for future work is to define a graphical
separation criterion along the lines of d-separation and t-separation
(Sullivant et al., 2010), which, given a mixed graph and a graphical
constraint, allows us to decide whether that constraint holds in the graph’s
model. Such a criterion may build on the restricted trek separation criterion
of Drton et al. (2020): our graphical constraints can be seen as recursively
nested determinants that obey additional structural restrictions.

Acknowledgments

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 883818).

## References

  * Cox et al. (2015) D. A. Cox, J. Little, and D. O’Shea.  _Ideals, Varieties, and Algorithms_.  Springer, New York, fourth edition, 2015.
  * Drton (2018) M. Drton.  Algebraic problems in structural equation modeling.  In _Advanced Studies in Pure Mathematics_ , pages 35–86. Mathematical Society of Japan, 2018.
  * Drton et al. (2020) M. Drton, E. Robeva, and L. Weihs.  Nested covariance determinants and restricted trek separation in Gaussian graphical models.  _Bernoulli_ , 26(4):2503–2540, 2020.
  * Foygel et al. (2012) R. Foygel, J. Draisma, and M. Drton.  Half-trek criterion for generic identifiability of linear structural equation models.  _The Annals of Statistics_ , 40(3):1682–1713, 2012.
  * Leung and Drton (2018) D. Leung and M. Drton.  Algebraic tests of general Gaussian latent tree models.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems 31 (NeurIPS 2018)_ , pages 6301–6310, 2018.
  * Robins (1986) J. Robins.  A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect.  _The Annals of Statistics_ , 7:1393–1512, 1986.
  * Shpitser et al. (2014) I. Shpitser, R. J. Evans, T. S. Richardson, and J. M. Robins.  Introduction to nested Markov models.  _Behaviormetrika_ , 41(1):3–39, 2014.
  * Shpitser et al. (2018) I. Shpitser, R. J. Evans, and T. S. Richardson.  Acyclic linear SEMs obey the nested Markov property.  In _Proceedings of the 34th Annual Conference on Uncertainty in Artificial Intelligence (UAI 2018)_ , 2018.
  * Spirtes et al. (2000) P. Spirtes, C. Glymour, and R. Scheines.  _Causation, Prediction, and Search_.  MIT Press, Cambridge, MA, second edition, 2000.
  * Sullivant et al. (2010) S. Sullivant, K. Talaska, and J. Draisma.  Trek separation for Gaussian graphical models.  _The Annals of Statistics_ , 38(3):1665–1685, 2010.
  * Tian and Pearl (2002) J. Tian and J. Pearl.  On the testable implications of causal models with hidden variables.  In _Proceedings of the 18th Conference on Uncertainty in Artificial Intelligence (UAI 2002)_ , 2002.
  * van Ommen and Mooij (2017) T. van Ommen and J. M. Mooij.  Algebraic equivalence of linear structural equation models.  In _Proceedings of the 33rd Annual Conference on Uncertainty in Artificial Intelligence (UAI 2017)_ , 2017.
  * Verma and Pearl (1991) T. S. Verma and J. Pearl.  Equivalence and synthesis of causal models.  In _Proceedings of the Sixth Conference on Uncertainty in Artificial Intelligence (UAI 1991)_ , 1991.
  * Weihs et al. (2017) L. Weihs, B. Robinson, E. Dufresne, J. Kenkel, K. Kubjas, M. Reginald, N. Nguyen, E. Robeva, and M. Drton.  Determinantal generalizations of instrumental variables.  _Journal of Causal Inference_ , 6, 12 2017.  doi: 10.1515/jci-2017-0009.