# A quintessence dynamical dark energy model from ratio gravity

Jackie C.H. Liu<EMAIL_ADDRESS>Institute for Advanced Study, Hong Kong
University of Science and Technology, Hong Kong  Department of Physics, Hong
Kong University of Science and Technology, Hong Kong

###### Abstract

Based on the work of ratio gravity developed in 2018, which postulates the
deformation of the cross ratio to associate with the physical model of
gravity, we develop a mechanism to generate dynamical dark energy - a
quintessence field coupled with gravity. Such model causes the dark energy
behaving differently in early and late time universe. In the radiation-
dominated-era and matter-dominated-era, the related analytical solutions of
the quintessence field have an interesting property - starting as a constant
field, then oscillating as the universe expands. By Markov Chain Monte Carlo
search of the parameter space with the local measurement (Type Ia supernovae)
in the Bayesian framework, the probed range of $H_{0}$ (within $1\sigma$)
overlaps the $H_{0}$ value inferred from Planck CMB dataset by $\Lambda$CDM
model.

Gravity, Cosmology

##  I INTRODUCTION

The accelerated expansion of the universe was discovered in 1998 by the
observation of supernova Riess  _et al._ (1998). The puzzle of dark energy is
one of the greatest problems in cosmology. Theorists propose different
explanations; for example, a simple one is the cosmological constant,
$\Lambda$, of the $\Lambda$CDM model, which is a widely accepted model because
it is part of the theoretical framework of general relativity, GR, and
consistent with the observations in great detail Ade  _et al._ (2016); Aghanim
_et al._ (2018). Alternative theories such as scalar field model of dark
energy Brax  _et al._ (2016) are also compelling cosmological models.

The measurements of the early universe and the local measurements (i.e. late-
time universe) Aghanim  _et al._ (2018); Riess  _et al._ (2016) probe the
Hubble constant for different values with increasing accuracy so both results
seem contradict, which is called the Hubble tension. It leads to many active
research because of the potential implication for new Physics for our
understanding of the gravity and cosmology Basilakos and Nesseris (2016).

There are several possible explanations to the tension such as statistical
fluke, or the emerging spatial curvature effect from the cosmological model of
the relativistic and nonlinear effect Bolejko (2018). The study of the Planck
data Xia  _et al._ (2013) shows that, one of possible solutions to the Hubble
tension, is that the equation of state of parameter of the dark energy, w, is
not equal to -1. Another proposed resolution is the early dark energy model,
EDE. EDE models drive the expansion of the early universe (usually in
radiation-dominated-era and/or matter-dominated-era), so that the expansions
of the universe in the early era and late time era behave differently. In the
work of Refs. Poulin  _et al._ (2018, 2019), one of the proposed axion models
initially freezes at constant field value, then evolves to oscillate after the
critical redshift so the EDE effect drives the early and late time of the
expansion periods differently. Poulin et al. Poulin  _et al._ (2019) analyse
the models against the Planck CMB dataset and Type Ia supernovae dataset to
show that the Hubble constant, $H_{0}$, probed by the model is consistent with
the $H_{0}$ inferred by Planck CMB dataset. Furthermore, another proposed dark
energy model is the acoustic dark energy model Lin  _et al._ (2019). It is
related to the dark fluid that a scalar field converts its potential energy to
kinetic energy during the matter-radiation equality.

The theory of ratio gravity, RG theory, is a newly developed theory Liu and
Wang (2018) that postulates the deformation of the cross ratio to associate
with the physical model of gravity in the framework of Newman-Penrose
formalism Newman and Penrose (1962). In the present work, we deploy a
different approach from Ref.Liu and Wang (2018) that we derive the physical
models of fermion and scalar fields from the core equations of the RG theory
in sections 2 and 3. We develop a generic framework to obtain the scalar
field, which has the property of the symmetry breaking for vacuum expected
value that is similar to the ordinary $\phi^{4}$ theory (in section 3).

In this paper, we explore the interesting property of a quintessence dynamical
dark energy model originated by the work of ratio gravity. We study the
prediction of the model in the the second part of the paper accordingly.

In section 4, by considering the scalar field as the quintessence dark energy
model with CDM, we show the correspondence to the $\Lambda$CDM model in
radiation-dominated and matter-dominated eras. In such eras, we found the
related analytical solutions of the field. The solutions have an interesting
property - starting as a constant field, then oscillating as the universe
expands. The models of axion dark energy Poulin  _et al._ (2018, 2019) suggest
similar scenario: the axion starts with ”frozen” phase then transits to the
oscillating phase.

In the last section, we perform the Markov Chain Monte Carlo search for the
parameters of the qCDM model - quintessence and CDM model - with the Pantheon
dataset - 1048 Type Ia supernovae (SNe Ia) Scolnic  _et al._ (2018). The
probed Hubble constant is approximately equals to 67$\pm$4 km/s/Mpc. Although
we probe the model parameters by the dataset of late time universe, the probed
range of $H_{0}$ (within $1\sigma$) surprisingly overlaps the $H_{0}$ deduced
from $\Lambda$CDM by Planck CMB dataset Aghanim  _et al._ (2018). Due to the
limited data analysis in this work, we make no conclusion to the possibility
for resolving Hubble tension by this model. Additional data analysis with more
dataset such as BAO is recommended in the future work.

##  II Introduction to the Framework

In this section, we introduce the principle of ratio gravity from the previous
work Liu and Wang (2018) and two core equations used throughout this paper. We
develop a new framework that relies on the basic principle of Ref. Liu and
Wang (2018), while modify the interpretation of the connection to gravity. We
explain the difference at the end of this section.

The definition of the cross ratio over Riemann sphere is:

$(z_{1},z_{2};z_{3},z)=\frac{(z_{3}-z_{1})(z-z_{2})}{(z_{3}-z_{2})(z-z_{1})},$

where $z_{1},z_{2},z_{3}$ are complex numbers of the poles, and $z$ is the
reference point over the Riemann sphere. A cross ratio consists of many
equivalent representations while represents the same value. The arbitrariness
of the same cross ratio allows four degrees of freedom because only four free
parameters for three movable poles. By Ref. Yoshida (1997), one of the
representations of the cross ratio is the hypergeometric differential equation
of three regular singular poles.

Since the hypergeometric differential equation can be written as a second
order linear differential equation in two-by-two matrix form and one can
express it as the integrable system Cassidy and Singer (2005):

$\partial_{\mu}Y=B_{\mu}Y.$

We can further introduce the gauge transformation to $Y$ and B matrices with
re-definitions of B matrices to yield:

$D_{\mu}Y=iB_{\mu}Y,$ (1)

where $D_{\mu}Y=\partial_{\mu}Y+i[\Lambda_{\mu},Y]$ and $\Lambda_{\mu}$ are
trace-less two-by-two Hermitian matrices. Eq. (1) is called $Y$ Equation. The
D operator obeys Leibniz rule for derivation. One can further transform $Y$
equation by tensoring a Hermitian map to yield the form:

$D_{ab}Y=iB_{ab}Y,$ (2)

Eq. (2) is the original form of $Y$ equation in Ref. Liu and Wang (2018) with
spinor index $ab$. Note that the Hermitian map can be associated with the
metric components according to Newman Penrose formalism (NP formalism Newman
and Penrose (1962)). The D operator can be defined in more general way as the
form $D_{\nu}=f^{\nu\mu}\partial_{\mu}$111$f^{\nu\mu}$ are complex functions.
as long as it is one-to-one corresponding to the D operator of $Y$ equation
(1) by associated automorphism. In this paper, we use capital index, e.g. $A$,
to denote the index abstractly to reserve the generalization for D operators.

The gauge transformation of Eq. (1) allows generating different
representations of the $Y$ equation trivially so it is an automorphism - the
transformed $Y$ equation and the original $Y$ equation are in the same space
(i.e. the same mathematical structure). Such automorphism allows we describe
the same cross ratio with different equivalent representations.

Galois transformation is introduced to provide another form of transformation
that respects the automorphism Liu and Wang (2018). The Galois transformation
is defined by a Galois operator, $\hat{\rho}$, that obeys:

$D\hat{\rho}=\hat{\rho}D,$ (3)

Eq. (3) is called the Galois equation. The definition of the Galois
transformation and Galois equation originate from the Galois differential
theory222Galois map $\pi:K\to K$, where $K$ is called differential field
extension. Cassidy and Singer (2005) \- the theory studies the Galois groups
of the differential equations. In the context of ratio gravity, we focus on
how Galois equation provides the transformation of $Y$ equation and the
related automorphism, so it requires no intensive knowledge of Galois
differential theory.

In section 3, we derive the equation of motion from $Y$ equation and introduce
a generic Galois operator that solves Galois equation and leads to the related
scalar field equation.

Unlike the previous work Liu and Wang (2018), we use original framework of
General Relativity instead of NP formalism. In previous work, the NP formalism
is connected to set of Galois equations via the introduction of Bianchi
constraints.

In this new framework (section 3), we first find the equation of motion of $Y$
equation (1), and the related scalar field equation(s) from Galois equation
(3) to define the associated Lagrangian of matter - $\mathcal{L}_{m}$. Then,
as the ordinary treatment of general relativity, we consider $\mathcal{L}_{m}$
as the source of gravity to define gravitational energy momentum tensor Wald
(1984), i.e. using Einstein equation as the constraint equation to fix the
degree of freedom of the metric.

In the context of RG, we use both Galois transformation and the continuous
transformation of the metric elements by general relativity (GR) to find the
space of transformed $Y$ equations, i.e. associated cross ratio
representations. Note that the automorphism in the context of RG is not the
same as the one in GR context - the general covariant transformations defined
as the automorphisms of fibre bundles; RG requires the automorphisms applying
to the space of $Y$ equations, i.e. the cross ratio representations.

##  III $Y$-Fermions and The Vacuum

In order to find the equation of motion of $Y$ equation (1) and the associated
Galois equation (3) in this section, we apply the gauge transformation to Eq.
(1), make use of Dirac equation, and find the related Galois equations. In the
middle and last parts of this section, we explain the interpretation of the
equation of motion of $Y$ equation and related Galois equation in the context
of quantum field theory. The purpose of this section is to define the physical
models in the RG context under the framework of Lagrangian.

By applying a gauge transformation to Eq. (1), one of the four components of
$Y$ matrix can be gauged out because there are 3 degrees of freedom of the
gauge in SU2. Therefore, there are four possible cases to choose the zeros of
the components of $Y$ matrix. We classify them as four categories: CAT 1$\to$4
of $Y$ matrix as follow

$\left(\begin{array}[]{cc}*&0\\\ &*\\\
\end{array}\right),\left(\begin{array}[]{cc}*&*\\\ 0&*\\\
\end{array}\right),\left(\begin{array}[]{cc}*&*\\\ &0\\\
\end{array}\right),\left(\begin{array}[]{cc}0&*\\\ &*\\\ \end{array}\right).$

In order to define the equation of motion by the eigen solutions of Eq. (1),
we introduce the parameterization to B and $\Lambda$ matrices:

$\displaystyle B_{\mu}$
$\displaystyle=\overset{1}{\Phi}\,p_{1\mu}\hat{e}+\overset{2}{\Phi}\,p_{2\mu}\hat{f}+\overset{3}{\Phi}\,p_{3\mu}\hat{h}+\overset{4}{\Phi}\,p_{4\mu}1_{2},$
(4) $\displaystyle\Lambda_{a\mu}$
$\displaystyle=\overset{5}{\Phi}\,p_{5\mu}\hat{e}+\overset{6}{\Phi}\,p_{6\mu}\hat{f}+\overset{7}{\Phi}\,p_{7\mu}\hat{h},$
(5)

where we use matrix structure of $sl_{2}$ algebras
($\hat{e},\hat{f},\hat{h}$), $p_{a\mu}$ are dimensionless parameters, and
$\overset{a}{\Phi}$ are complex functions. We define $y$ for the column matrix
of three dimensions to represent the non-zero components of $Y$ matrix. We can
re-write the $Y$ equation for CAT 1 $\to$ CAT 4 as the following form:

$i\partial_{\mu}y=P_{\mu}y,$ (6)

where $P_{\mu}$ are the three-by-three matrices. There are constraints of
{$\overset{a}{\Phi}$} needed to be satisfied to obtain Eq. (6); for instance,
the explicit form of $P_{\mu}$ matrices for CAT 1 is

$\small
P_{\mu}=\left(\begin{array}[]{ccc}-\overset{4}{\Phi}p_{4\mu}&0&-\overset{3}{\Phi}p_{3\mu}\\\
-\overset{2}{\Phi}p_{2\mu}&\overset{3}{\Phi}p_{3\mu}-\overset{4}{\Phi}p_{4\mu}-2\overset{7}{\Phi}p_{7\mu}&-\overset{2}{\Phi}p_{2\mu}\\\
-\overset{3}{\Phi}p_{3\mu}&0&-\overset{4}{\Phi}p_{4\mu}\\\
\end{array}\right),$ (7)

with $\overset{1}{\Phi}=\overset{5}{\Phi}=\overset{6}{\Phi}=0$.

The eigen-values of $P_{\mu}$ matrices of Eq. (6) for CAT 1,2,3,4 are
{$\pm\phi_{1}p_{\mu},\phi_{2}p_{\mu}$},
{$\mp\phi_{1}p_{\mu},-\phi_{2}p_{\mu}$},
{$\pm\phi_{2}p_{\mu},\phi_{1}p_{\mu}$}, and
{$\mp\phi_{2}p_{\mu},-\phi_{1}p_{\mu}$} respectively, where $p_{\mu}$ is the
dimensionless momentum constructed by parameters $p_{a\mu}$, and $\phi_{1}$
and $\phi_{2}$ are the linear combinations of {$\overset{a}{\Phi}$}. Because
of the first order differential operator of Eq. (6), we make use of Dirac
equation in momentum space {$L(p),R(p)$} to obtain the equation of motion of
$y$: 333The Dirac-slash-notation operators $\not\\!p$ and $\not\\!\partial$
contain the projector matrix for $l\leftrightarrow r$ implicitly.
$\not\\!p\phi^{*}L=m_{y}\phi R$ and $\not\\!p\phi_{2}R=m_{y}\phi^{\dagger}L$
are the constraints satisfied at real classical expected value of $\phi$.

$\displaystyle i\not\\!\partial L=\not\\!p\,\phi^{*}L=m_{y}\,\phi R,$
$\displaystyle i\not\\!\partial R=\not\\!p\,\phi_{2}R=m_{y}\,\phi^{\dagger}L,$
(8)

where $m_{y}$ is the mass coupling, $L=L(p)l,R=R(p)r$ and $l^{*}$ are the
doublet of first and second eigen-vectors, and $r$ is the third eigen-vector
of $y$ respectively444The definition of $l^{*},r$ to correspond to eigen-
vectors of $y$ is only conventional because of the standard model framework;
one can employ $l,r^{*}$ with the re-definition of $\phi\to\phi^{*}$., and we
denote $\phi$ as the doublet form of {$\phi_{1},\phi_{2}$}. We naturally
define the equation of motion of $y$ as the equation of motion for $Y$
fermion. The equation of motion cannot be solved because the value of $\phi$
are not constrained, so we rely on Galois equation Eq. (3) to fix it next.

In the context of quantum field, we interpret that the excitation of multiple
$Y$ fermions by the Dirac Lagrangian associated with Eq. (8) is corresponding
to the set of multiple representations of the related $Y$ equation. It is
merely the interpretation to relate the context of quantum field from the RG
theory’s perspective.

Cassidy Cassidy and Singer (2005) defines the Galois map (automorphism $\pi$)
to transform as $\pi:x\to y$, where $x$ and $y$ are elements of the space
constructed by $Y$ matrix and the derivatives of $Y$ matrix, and
$\pi\partial=\partial\pi$ as Eq. (3). We define a generic Galois operator,
$\hat{\rho}$, similarly: $\hat{\rho}:=X^{A}D_{A}$, where $\hat{\rho}$
satisfies Galois equation, i.e.

$D_{C}\,\hat{\rho}\,Y=\hat{\rho}\,D_{C}\,Y,$ (9)

and $X^{A}$ are two-by-two-matrix functions because the operator D acts on
two-by-two matrices. Certainly, one can define a more complicated Galois
operator (e.g. higher derivative operator) with the potential cost of less
solvability of Galois equation. To ensure $\hat{\rho}$ being associated with
automorphism, we use the exponential map: $exp(\epsilon\hat{\rho})$, such that
$Y$ equation (1) transforms invariantly and infinitesimally by $\epsilon$ if
Eq. (9) is satisfied.

Because of the parameterization for $Y$ equation, the Galois equation (9) can
be generally expressed as (for CAT 1 and 4):

$\displaystyle\omega^{A}\left(h_{\mathcal{C}}\partial_{A}\left(\overset{3}{\Phi}\right)-h_{A}\partial_{\mathcal{C}}\left(\overset{3}{\Phi}\right)\right)-\overset{3}{\Phi}h_{A}\partial_{\mathcal{C}}\left(\omega^{A}\right)=0,$
$\displaystyle\omega^{A}f_{\mathcal{C}}\partial_{A}\left(\overset{2}{\Phi}\right)-f_{A}\omega^{A}\partial_{\mathcal{C}}\left(\overset{2}{\Phi}\right)+$
$\displaystyle\overset{2}{\Phi}\left(2\overset{7}{\Phi}\omega^{A}\left(f_{A}j_{\mathcal{C}}-j_{A}f_{\mathcal{C}}\right)-f_{A}\partial_{\mathcal{C}}\left(\omega^{A}\right)+2\overset{3}{\Phi}\omega^{A}\left(h_{A}f_{\mathcal{C}}-f_{A}h_{\mathcal{C}}\right)\right)+$
$\displaystyle\mathit{f}^{A}\left(2\overset{3}{\Phi}\overset{7}{\Phi}h_{A}j_{\mathcal{C}}-h_{A}\partial_{\mathcal{C}}\left(\overset{3}{\Phi}\right)+h_{\mathcal{C}}\partial_{A}\left(\overset{3}{\Phi}\right)\right)-\overset{3}{\Phi}h_{A}\partial_{\mathcal{C}}\left(\mathit{f}^{A}\right)=0,$
(10)

where $j_{A},h_{A},f_{A}$ denote the parameters $ip_{7A},ip_{3A},ip_{2A}$
respectively, $\overset{4}{\Phi}$ is zero, and
$X(x)^{A}=\omega(x)^{A}1_{2}+\mathit{f}(x)^{A}\hat{f}$. For CAT 2 and 3,
equations (10) are the same form by the transformation:
$\mathit{f}(x)^{A}\to\mathit{e}(x)^{A}$ and $f_{A}\to e_{A}$. We notice the
equation above can be realized in a symbolic form as:

$\partial\Phi=(\Phi)+(\Phi\Phi),$ (11)

where ($\Phi$) and ($\Phi\Phi$) denote the terms with the coefficients for the
powers of ($\overset{a}{\Phi}$) and ($\overset{a}{\Phi}\overset{b}{\Phi}$)
respectively, the equation is likely in the form for scalar field(s) with non-
zero vacuum expected value, vev. The rest of this section is to prove this
observation and construct the associated symmetry-breaking Lagrangian.

By applying the rest frame condition onto Eq. (10) for the dimensionless
momentum, $p=(p_{O},0)$, of $y$ fermion, we obtain:

$\displaystyle\phi_{1}\partial_{\mathcal{C}}\left(\mathit{f}^{O}\right)+\mathit{f}^{O}\partial_{\mathcal{C}}\left(\phi_{1}\right)=0,$
$\displaystyle\mathit{f}^{O}\left(\mu\phi_{1}^{2}+i\partial_{O}\left(\phi_{1}\right)\right)+i\phi_{1}\partial_{O}\left(\mathit{f}^{O}\right)-i\mathit{f}\cdot{\nabla}\left(\phi_{1}\right)=0$
(12)

where $O$ denotes the time-axis-index, $O\neq C$, $\mathit{f}\cdot\nabla$
denotes directional derivate $\mathit{f}^{A}\partial_{A}$, and $\mu:=p_{O}$ is
the dimensionless constant of the theory. Eq. (12) is only a specific solution
of Eq. (10) when we consider the case of singlet $\phi_{1}$ solution. (The
doublet equation is not covered in this paper.) In the rest of this paper, we
denote $\phi$ as the singlet field.

In order to apply to a specific coordinate system for cosmology, we consider
the $D_{O}$ operator by $D_{O}=a(t)\partial_{t}$ for FRW cosmology. The
Laplacian of the time-dependent-only $\phi$ is555We use the degree of freedom
of $\mathit{f}^{C}$, $C\neq O$, to remove the directional derivate term.

$\ddot{\phi}=-\frac{i\mu\phi^{2}\dot{a}}{a^{2}}-\frac{2\mu^{2}\phi^{3}}{a^{2}}-\frac{3i\mu\phi^{2}\dot{\theta}}{a\theta}+\frac{2\phi\dot{\theta}^{2}}{\theta^{2}}-\frac{\phi\ddot{\theta}}{\theta},$
(13)

where $\theta$ denotes $\mathit{f}^{O}$, and the associated Lagrangian is

$\displaystyle-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\chi_{1}\partial_{\nu}\chi_{1}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\chi_{2}\partial_{\nu}\chi_{2}-$
$\displaystyle(i\chi_{1}^{3}-3\chi_{2}\chi_{1}^{2}-3i\chi_{2}^{2}\chi_{1}+\chi_{2}^{3})C+$
$\displaystyle\chi_{1}^{2}A-\chi_{2}^{2}A+\chi_{1}\chi_{2}B-\frac{\mu^{2}\left(\chi_{1}^{4}+\chi_{2}^{4}\right)}{2a^{2}}+$
$\displaystyle\chi_{1}^{2}\left(\frac{3\mu^{2}\chi_{2}^{2}}{a^{2}}+\frac{\mu\chi_{2}\dot{a}}{a^{2}}\right)-\frac{\mu\chi_{2}^{3}\dot{a}}{3a^{2}},$
(14)

where $\phi$ is expressed in real and imaginary parts:
$\phi=\chi_{1}-i\,\chi_{2}$, $B=-2iA$, $C=\frac{\mu\dot{\theta}}{2a\theta}$
and $A=\frac{\dot{\theta}^{2}}{2\theta^{2}}-\frac{\ddot{\theta}}{4\theta}$.
The Lagrangian above is not yet the physical model we look for. The
problematic complexness does not respect the Hermiticity of Lagrangian. So, we
add the Hermitian conjudge terms. Because the Hermitian map of D operator (8),
the associated $Y$ fermion respects the symmetry of positive and negative
energies. Therefore, we model the Lagrangian terms associated with the Galois
equation in the current theory as
$\mathcal{L}_{\phi}=\mathcal{L}_{+\phi}+\mathcal{L}_{-\phi}$ such that it
respects positive-negative-vev-symmetry, just like the ordinary $\phi^{4}$
theory respects $Z_{2}$ symmetry. Obviously, this artificial symmetry breaks
down if the $Y$-fermionic sector does not obey such symmetry. We consider this
possibility to be the future development. Finally, we obtain the effective
potential $V_{\chi}$:

$\frac{\mu^{2}\left(\left(\chi^{*}\right)^{4}+\chi^{4}\right)}{4a^{2}}-\frac{1}{2}\left(\chi^{2}+\left(\chi^{*}\right)^{2}\right)m_{\chi}(\theta)^{2},$
(15)

where $m_{\chi}(\theta)^{2}:=A$.

In this section, we show how to obtain the fermionic model of theory, i.e. $Y$
fermion, associated with $Y$ equation (1) for CAT 1$\to$4, and the symmetry-
breaking scalar field potential (15) by the generic Galois operator and Galois
equation.

##  IV The Quintessence Field and Cosmological Model

In this section, we apply the framework of previous section to the application
of cosmology - simplify the $\chi$ potential (15) to construct the
quintessence field model that coupled to gravity. We show that such qCDM model
corresponds to well-accepted $\Lambda$CDM with a derivation mechanism. Brief
comparison to several established models Poulin  _et al._ (2018, 2019) is
covered.

The vacuum expected value of $\chi$ potential (15) is
$\frac{m_{\chi}(\theta)}{\mu/a}$. It is not fixed because of the degree of
freedom by $m_{\chi}(\theta)$, and it is dynamical as the scale factor varies.
By requiring the vev of $\chi$ fixed, the $m_{\chi}(\theta)$ term should be
proportional to $1/a$ so the vev becomes $\frac{m_{\chi}}{\mu}$, and the
$\chi$ potential becomes

$\frac{\mu^{2}\left(\left(\chi^{*}\right)^{4}+\chi^{4}\right)}{4a^{2}}-\frac{1}{2a^{2}}\left(\chi^{2}+\left(\chi^{*}\right)^{2}\right)m_{\chi}^{2},$
(16)

where $m_{\chi}$ is a constant parameter. In the rest of the paper, we
consider the simple case to define the quintessence field, that $\chi$ is a
real scalar field which recovers as $\phi^{4}$ potential at $a=1$, and it has
the minimum degree of freedom needed to solve Friedmann equations.

The quintessence field potential is re-written as

$\frac{\mu^{2}\chi^{4}}{4a^{2}}-\frac{m_{\chi}^{2}\chi^{2}}{2a^{2}}.$ (17)

We consider the Lagrangian of the quintessence666$Z_{\chi}$ is the overall
coupling strength between gravity and $\chi$ field.,
$\mathcal{L}_{q}:=Z_{\chi}\mathcal{L}_{\chi}$, as the source of gravity to
define gravitational energy momentum tensor for the quintessence field, and
apply the usual variation on $\mathcal{L}_{q}$ to get the related density and
pressure of the quintessence field777Prime denotes the derivative with respect
to $a^{2}$. Wald (1984):

$\displaystyle\rho_{q}$
$\displaystyle=\frac{\mu^{2}q^{4}}{4a^{2}}-\frac{m_{\chi}^{4}}{4a^{2}\mu^{2}}+\frac{\mu
q^{3}m_{\chi}}{a^{2}}+\frac{q^{2}m_{\chi}^{2}}{a^{2}}+\frac{1}{2}Z_{\chi}\dot{q}^{2},$
$\displaystyle p_{q}$
$\displaystyle=\frac{\mu^{2}q^{4}}{4a^{2}}-2\mu^{2}q^{3}q^{\prime}-\frac{m_{\chi}^{4}}{4a^{2}\mu^{2}}+\frac{\mu
q^{3}m_{\chi}}{a^{2}}+\frac{q^{2}m_{\chi}^{2}}{a^{2}}-$ $\displaystyle 6\mu
q^{2}m_{\chi}q^{\prime}-4qm_{\chi}^{2}q^{\prime}+\frac{1}{2}Z_{\chi}\dot{q}^{2},$

where we expand the $\chi$ field around vev, $\chi=\frac{m_{\chi}}{\mu}+q$, so
we can deploy weak field limit next. We notice that $V_{\chi}(a^{2})$
contributes because of the variation, and absorb $Z_{\chi}$ factor for terms,
$Z_{\chi}\mu^{2}\to\mu^{2}$ and $Z_{\chi}m_{\chi}^{2}\to m_{\chi}^{2}$. With
the re-definition of the constants, the $Z_{\chi}$ factor is merely the
rescaling factor for the time/energy scale between the quintessence and $\chi$
fields. We further assume the validity of weak-field-limit, i.e. quadratic-
terms-dominated, to yield:

$\displaystyle\rho_{q}$
$\displaystyle=-\frac{m_{\chi}^{4}}{4a^{2}\mu^{2}}+\frac{q^{2}m_{\chi}^{2}}{a^{2}}+\frac{1}{2}Z_{\chi}\dot{q}^{2},$
$\displaystyle p_{q}$
$\displaystyle=-\frac{m_{\chi}^{4}}{4a^{2}\mu^{2}}+\frac{q^{2}m_{\chi}^{2}}{a^{2}}-4qm_{\chi}^{2}q^{\prime}+\frac{1}{2}Z_{\chi}\dot{q}^{2}.$
(18)

In order to study the dark energy behavior, we define the dark energy density
parameter, $\xi:=\Omega_{q}$, so

$\xi=-\frac{\alpha^{4}}{12a^{2}\lambda
H_{0}^{4}}+\frac{\alpha^{2}q^{2}}{3a^{2}H_{0}^{4}}+\frac{\lambda
Z_{\chi}\dot{q}^{2}}{6\mu^{2}H_{0}^{4}}$ (19)

where $\mu^{2}:=\frac{\lambda M_{p}^{2}}{H_{0}^{2}}$ and
$m_{\chi}:=\frac{\alpha M_{p}}{H_{0}}$ for the ease of parameter probing next.
By Friedmann equations (without curvature k and cosmological constant terms),
we have

$\displaystyle\dot{\xi}=$ $\displaystyle
H(\frac{\alpha^{6}}{\lambda}+\frac{2aFH_{0}^{2}\ddot{a}+FH_{0}^{4}a^{2-r}(4\xi
a^{r}+\Omega_{m}(3w+1))}{\mu^{2}q}-4\alpha^{4}q^{2})/$
$\displaystyle(6a^{2}\alpha^{2}H_{0}^{4}),$

where $F:=(2\alpha^{2}\mu^{2}q+a^{2}\lambda Z_{\chi}\ddot{q})$, $w$ is the
equation of state of the matter component, and $r=3,4$ for matter and
radiation-dominated eras respectively. We found if we define the equation of
motion for q

$\ddot{q}=-\frac{2\alpha^{2}\mu^{2}q}{\mathit{a}^{2}\lambda Z_{\chi}}$ (20)

then, when $F=0$, the weak field limit is valid up to a long period of
cosmological time span. The Friedmann equations together with Eq. (20) can be
expressed as

$\xi=\frac{H^{2}}{H_{0}^{2}}-\frac{\Omega_{m}}{a^{r}},\dot{\xi}=\frac{\alpha^{2}H\left(\alpha^{2}-4\lambda
q^{2}\right)}{6a^{2}\lambda
H_{0}^{4}},\ddot{q}=-\frac{2\alpha^{2}\mu^{2}q}{a^{2}\lambda Z_{\chi}}.$ (21)

Eq. (21) are the equations of the qCDM model of the quintessence theory.

The $\Lambda$CDM model is effective and supported by many observations. We
need to verify the validity of qCDM model analytically against $\Lambda$CDM
model. It is clear if $\dot{\xi}$ is zero, i.e. $\alpha^{2}=4\lambda q^{2}$,
then we recover the case of the cosmological constant and the first equation
of (21) is simply the Friedmann equation with the cosmological constant so the
$\Lambda$CDM correspondence is satisfied. Therefore, $q$ must be approximately
constant and equals to $\frac{\alpha}{2\sqrt{\lambda}}$ in order to justify
the validity of $\Lambda$CDM correspondence. It can be achieved by the
quintessence field staying approximately constant for long period of time or
oscillating very slowly.

In both radiation-dominated-era and matter-dominated-era, i.e. $Log(a)\propto
Log(t)$, we can solve $q$ analytically. In matter-dominated-era,

$q=(c_{2}+c_{1}\omega T)\sin(\omega T)+(c_{1}-c_{2}\omega T)\cos(\omega T),$
(22)

where $\omega=\frac{3\sqrt{2}\alpha\mu}{\sqrt{\lambda}\sqrt{Z_{\chi}}}$,
$T=\sqrt[3]{t/t_{0}}$, $c_{1},c_{2}$ are integration constants, and $t_{0}$ is
the present time, $a(t_{0})=1$; in radiation-dominated-era,

$q=c_{3}T\,J_{1}(\omega T)+c_{4}T\,Y_{1}(\omega T),$ (23)

where $J_{1}$ is a Bessel function of the first kind, $Y_{1}$ is a Bessel
function of the second kind,
$\omega=\frac{2\sqrt{2}\alpha\mu}{\sqrt{\lambda}\sqrt{Z_{\chi}}}$,
$T=\sqrt{t/t_{0}}$, and $c_{3},c_{4}$ are integration constants. Both
analytical expressions lead to the constant mode as $t\to 0$, so in the early
universe, the quintessence field is asymptotically constant and later evolves
to oscillate.

The Eq. (23) is similar to the scalar field model in Ref. Marsh and Ferreira
(2010) that showed ultra-light scalar fields affect growth of structure in the
Universe as well as the expansion rate. In the limit that, $a(t)\propto
t^{p}$, i.e. in both radiation-dominated and matter-dominated eras, the
analytical form of the axion-like particles is similar to (but not the same
as) Eq. (23).

Given that the quintessence is nearly constant, $q_{c}$ , $\xi$ is solved

$\xi=\xi_{0}-\frac{C}{a^{2}},$ (24)

where $\xi_{0}$ is an effective cosmological constant term and
$C=\frac{\alpha^{4}-4\alpha^{2}\lambda q_{c}^{2}}{12H_{0}^{4}\lambda}$ as a
constant parameter associated with the term scaling as the spatial curvature
(i.e. $a^{-2}$). We have shown by the analytical form of the quintessence
field that the near-constant-approximation is applicable in the early universe
(in radiation-dominated-era and matter-dominated-era), and then the field
oscillates shown in FIG. 1; the dark energy density parameter $\xi$ has an
effective cosmological constant term in Eq. (24).

Figure 1:  Sample plot of the analytic formula of quintessence in the matter-
dominated-era: $\alpha=2,\lambda=1,Z_{\chi}=1$. $t_{0}=1$ is the present time.

Interestingly, the models of dark energy Poulin  _et al._ (2018, 2019) suggest
similar scenario (but not exactly the same): the axion starts with ”frozen”
phase as the cosmological constant then transits to the oscillating phase,
FIG.1 of Ref. Poulin  _et al._ (2018).

##  V Data Analysis

In our data analysis section, we first identify a smaller set of parameter
spaces of Eq. (21) that reduces the qCDM model to an effective and simplified
qCDM version; following the same procedure as in Poulin  _et al._ (2019); Wang
(2018), we use the Pantheon dataset of 1048 SNe Ia Scolnic  _et al._ (2018) to
probe the model parameters with $\Omega_{m}$ dominated by matter component
only.

The original parameters of qCDM are
{$H_{0},\Omega_{m},Z_{\chi},\lambda,\alpha,q_{0},q_{1}$}, where $q_{0}$ and
$q_{1}$ denote the $q(t_{0})$ and $\dot{q}(t_{0})$ of the present day. We fix
$Z_{\chi}=\mu^{2}$ because it allows that $\lambda$ is in scale within unity
and $\alpha$ is in scale of $H_{0}$. We also assume the effectiveness of
$\Lambda$CDM in the present day such that $\dot{\xi}(t_{0})$ and $q_{1}$ are
effectively zero so $4q_{0}^{2}=\alpha^{2}\lambda$. The parameters of our
simplified qCDM model are {$H_{0},\Omega_{m},\lambda,\alpha$}.888The base unit
of $\alpha$ and q is set as 73.9 km/s/Mpc.

The sampling by the numerical solving for Eqs. (21) consumes the computation
resources seriously; therefore, we impose the prior-assumption. By applying
Markov Chain Monte Carlo, MCMC, parameters searching999An open sourced MCMC
Mathematica package with modification for the our model
(https://github.com/joshburkart/mathematica-mcmc). and assuming the flat
priors101010We use the preferred prior from the known cosmological parameter:
$0.24<\Omega_{m}<0.36$, $57<H_{0}<85$ km/s/Mpc. on
{$H_{0},\Omega_{m},\lambda,\alpha$}, we perform the initial MCMC exploration
on binned data (40 data-points), then we identify the preferred prior region
for $\alpha\leq 0.25$, whereas the small-valued $\alpha$ is theoretically
suggested because $\xi$ becomes effectively the cosmological constant as
$\alpha\to 0$ by Eq. (21). The full MCMC run on the Pantheon dataset (1048 SNe
Ia) Scolnic  _et al._ (2018) shows the region of convergence from the 1D and
2D posterior distributions in FIG. 2 with Gelman-Rubin criterion $R-1<0.024$.

Figure 2:  1D and 2D posterior distributions of $H_{0},\alpha,\Omega_{m}$
probed by the late time measurement (Pantheon dataset of 1048 SNe Ia). The
mean $H_{0}$ ($\pm 1\sigma$) of qCDM model is approximately equals to 67$\pm$4
km/s/Mpc; the referenced $H_{0}$ of Planck Aghanim  _et al._ (2018) and the
referenced $H_{0}$ of the late time measurements Riess  _et al._ (2016) are
shown in red and blue respectively.

The $\chi^{2}_{min}$ in the MCMC run is $1035.7$, that is slightly better than
$1036.5$ from the $\chi^{2}_{min}$ of $\Lambda$CDM probe Wang (2018). The mean
$H_{0}$ ($\pm 1\sigma$) of qCDM model is 67$\pm$4 km/s/Mpc, that overlaps the
estimated range (67.4$\pm$0.5 km/s/Mpc) by Planck measurement Aghanim  _et
al._ (2018). We found the acceptable range of parameter $\lambda$ is board; we
report the mean-best-fit-parameters111111Taking the mean values of parameters
of the top 2% best-fit parameters within the 1$\sigma$ region of $H_{0}$,
$\Omega_{m}$, and $\alpha$. as {$H_{0}=$67.2
km/s/Mpc$,\Omega_{m}=0.285,\lambda=0.698,\alpha=0.132$}.

Finally, we use the Eq. (24) and the mean-best-fit-parameters to check briefly
if the $C$ term of qCDM model is consistent with the CMB power spectrum.
Without a complete probe of qCDM parameters, we only change the $\Lambda$
term, cosmological constant of $\Lambda$CDM, by the Eq. (24). The value of
$q_{c}$ is $0.071$, and the order of magnitude of $C$ is $-4.9$. We obtain the
values of $\chi^{2}$ of CMB power spectrum against the dataset of PlanckTT and
WMAPTT121212PlanckTT dataset is corresponding to 111 data-points of $2\leq
l\leq 2000$, and WMAPTT dataset is corresponding to 45 data-points of $2\leq
l\leq 1150$ from the open-source-CMB tool
(http://www2.iap.fr/users/pitrou/cmbquick.htm). The $\chi^{2}$ of PlanckTT and
WMAPTT are only shifted by 0.02 and 0.07 respectively ($\chi^{2}$ for PlanckTT
and WMAPTT are 493.37, 42.54 respectively). However, the complete parameter
probe of qCDM against CMB power spectrum is not covered in this work.

##  VI Discussion

In this work, we introduce the framework of ratio gravity that postulates the
transformation of cross ratio is related to different representations of the
associated fermion and scalar models. The theory provides the mechanism to
generate the symmetry-breaking scalar fields naturally, that leads to the
quintessence field to drive the dark energy behaving dynamically.

The presented qCDM model can reproduce the $\Lambda$CDM model with a
derivation mechanism. The data analysis of the model with the supernovae
dataset suggests $H_{0}$ = 67$\pm$4 km/s/Mpc, which is aligned with the latest
Planck observation Aghanim  _et al._ (2018). Yet, a further analysis with
complete set of qCDM parameters against CMB and BAO dataset is suggested.

The theory can be extended to the domain of complex singlet and doublet models
of the scalar field. As the mass scale of the quintessence field is in $H_{0}$
suggested in the section 5, the possibility of light-massive boson because of
the complex phase of the singlet model is worth to be studied.

###### Acknowledgements.

JCHL would like to thank Professor Wang Yi and Professor Henry Tye for the
valuable comments and advices. This research was conducted using computational
resources at the Institute for Advanced Study, Hong Kong University of Science
and Technology.

## References

  * Riess  _et al._ (1998) A. G. Riess _et al._ , Astron.J. 116, 1009 (1998), astro-ph/9805201.
  * Ade  _et al._ (2016) P. Ade _et al._ , Astronomy and Astrophysics 594, A13 (2016).
  * Aghanim  _et al._ (2018) N. Aghanim _et al._ ,  (2018), arXiv:1807.06209v2[astro-ph].
  * Brax  _et al._ (2016) P. Brax, C. Burrage, C. Englert, and M. Spannowsky, Physical Review D 94 (2016), 10.1103/PhysRevD.94.084054.
  * Riess  _et al._ (2016) A. G. Riess _et al._ , The Astrophysical Journal 826 (2016), 10.3847/0004-637X/826/1/56.
  * Basilakos and Nesseris (2016) S. Basilakos and S. Nesseris, Physical Review D 94 (2016), 10.1103/PhysRevD.94.123525.
  * Bolejko (2018) K. Bolejko, Physical Review D 97 (2018), 10.1103/PhysRevD.97.103529.
  * Xia  _et al._ (2013) J.-Q. Xia, H. Li, and X. Zhang, Physical Review D 88 (2013), 10.1103/PhysRevD.88.063501.
  * Poulin  _et al._ (2018) V. Poulin, T. L. Smith, D. Grin, T. Karwal, and M. Kamionkowski, Physical Review D 98 (2018), 10.1103/PhysRevD.98.083525.
  * Poulin  _et al._ (2019) V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Physical Review Letters 122 (2019), 10.1103/PhysRevLett.122.221301.
  * Lin  _et al._ (2019) M.-X. Lin, G. Benevento, W. Hu, and M. Raveri, Physical Review D 100 (2019), 10.1103/PhysRevD.100.063542.
  * Liu and Wang (2018) J. C. Liu and Y. Wang, Phys.Rev.D 98 (2018).
  * Newman and Penrose (1962) Newman and R. Penrose, J.Math.Phys 3, 566 (1962).
  * Scolnic  _et al._ (2018) D. M. Scolnic _et al._ , The Astrophysical Journal 859 (2018).
  * Yoshida (1997) M. Yoshida,  _Hypergeometric Functions, My Love_ (Springer Fachmedien Wiesbader, Germany, 1997).
  * Cassidy and Singer (2005) Cassidy and Singer, arxiv.org (2005), arXiv:math/0502396[math.CA].
  * Wald (1984) R. M. Wald,  _General Relativity_ (The University of Chicago Press, 1984).
  * Marsh and Ferreira (2010) D. J. E. Marsh and P. G. Ferreira, Physical Review D 82 (2010), 10.1103/PhysRevD.82.103528.
  * Wang (2018) D. Wang, Physical Review D 97 (2018), 10.1103/PhysRevD.97.123507.


11institutetext: School of Mathematics, University of Leeds, Leeds LS2 9JT,
UK. Web: http://maths.leeds.ac.uk/~kisilv/

Molecular dynamics and particle methods Semiclassical theories and
applications Molecular dynamics and other numerical methods

# Comment on ‘Do we have a consistent non-adiabatic quantum-classical
mechanics?’

Vladimir V. Kisil On leave from Odessa University<EMAIL_ADDRESS>
###### Abstract

We argue with claims of the paper [1] that the quantum-classic bracket
introduced in [4] produces “artificial coupling” and has “genuinely classical
nature”.

###### pacs:

02.70.Ns

###### pacs:

03.65.Sq

###### pacs:

31.15.Qg

##  1 Introduction

This is a comment on the paper [1], which evaluates the quantum-classical (QC)
bracket:

$\displaystyle{}[K_{1},K_{2}]_{qc}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathrm{i}h}[K_{1},K_{2}]+\frac{1}{2}\left(\\{K_{1},K_{2}\\}-\\{K_{2},K_{1}\\}\right)$
(1)
$\displaystyle-\left.\mathrm{i}\partial_{h_{2}}[K_{1},K_{2}]\right|_{h_{2}=0},$

introduced in paper [4, (26)]. The authors in [1] claimed that the QC bracket
(1) exhibits:

  * •

an artificial coupling property (i.e., coupling between the subsystems in the
absence of an interaction);

  * •

a genuinely classical nature (i.e., the apparent mixed quantum classical form
reduces to a purely classical form for both subsystems).

The assessment in [1] oversaw the following points:

  1. 1.

QC bracket (1) is the image of the universal bracket [4, (22)]:

$\left\\{\\!\left[k_{1},k_{2}\right]\\!\right\\}=(k_{1}*k_{2}-k_{2}*k_{1})(\mathcal{A}_{1}+\mathcal{A}_{2}),$
(2)

under QC representation [4, (20)]. The universal bracket consists of
convolution commutator and antiderivative operators [4, (12)]. QC bracket
requires consideration of the first jet space [4]: the bracket is determined
not only by their values of observables at $h_{2}=0$ but also by the values of
their first derivative with respect to $h_{2}$ at zero (see the last term in
(1)).

  2. 2.

The derivation QC bracket (1) is independent of p-mechanisation procedure
introduced in [4, (23)]:

$\displaystyle q_{j}$ $\displaystyle\mapsto$ $\displaystyle
Q_{j}=\delta^{\prime}_{x_{j}}(g_{1};g_{2}),$ (3) $\displaystyle p_{1}$
$\displaystyle\mapsto$ $\displaystyle
P_{j}=\chi^{\prime}_{s_{k}}(s_{1}+s_{2})*\delta^{\prime}_{y_{j}}(g_{1};g_{2}),\quad$
(4)

where $j=1$, $2$ and $k=3-j$.

Here _p-mechanisation_ [3, § 3.3], as an analog of quantisation, is a
prescription how to build p-mechanical observables out of classical ones. It
may not be very explicit in [4], but the deduction of the bracket (1) is
compatible with different choices of p-mechanisation, however the value of the
bracket will be different, see Exs. 3 and 4 below. To illustrate this in the
present comment we use p-mechanisation given by the Weyl (symmetric) calculus
based on the following correspondence, cf. [4, (23)], [1, (19)] and (3)–(4):

$q_{j}\mapsto Q_{j}=\delta^{\prime}_{x_{j}}(g_{1};g_{2}),\quad p_{j}\mapsto
P_{j}=\delta^{\prime}_{y_{j}}(g_{1};g_{2}),\quad$ (5)

Then the quantum-quantum image of the universal bracket (2) of the respective
coordinate and momentum observables is:

$[Q_{j},P_{j}]_{qq}=\frac{h_{1}+h_{2}}{h_{k}}I,\qquad k=3-j.$ (6)

Now we review the above two claims from the paper [1].

##  2 Artificial coupling property

There is the following claim in [1, 3001-p3]: “It must be underlined that eq.
(16) describes an artificial interaction even if the two systems are not
coupled by the Hamiltonian.” This coupling property is attributed to the fact
that quantum-quantum bracket in [4, (25)] and [1, (16)] always contains both
Planck’s constants $h_{1}$ and $h_{2}$, which are generated by the presence of
both antiderivative operators in the definition of universal bracket (2).

######  Example 1.

In order to exam the claim let us consider an uncoupled Hamiltonian
$H(q_{1},p_{1},q_{2},p_{2})=H_{1}(q_{1},p_{1})+H_{2}(q_{2},p_{2})$. The
p-mechanisation (as well as quantisation) is a linear map [3, § 3.3], thus
this uncoupled structure will be preserved. Let $\hat{B}$ be an observable
depending only from $\hat{X}_{2}$ and $\hat{D}_{2}$, thus it will commute with
$H_{1}$. Therefore the commutator of $B$ and $H$ will be the same as $B$ and
$H_{2}$. The QC bracket is an image under a representation of the usual
commutator, thus the universal bracket (2) of $B$ and $H$ will be the same as
$B$ and $H_{2}$. Consequently the $\hat{H}_{1}$ will not affect the dynamics
of such an observable $\hat{B}$.

Therefore there is no coupling in the following meaning: arbitrary change of
the Hamiltonian $H_{1}$ of the first subsystem will not affect dynamics of any
observable build from coordinates and momenta of the second system only.

##  3 Genuinely classical nature

The paper [1, 3001-p3] said “In ref. [8] it was suggested that the dynamical
equation (16), in the limit $h_{1}=h$ and $h_{2}\rightarrow 0$, yields a QC
dynamics.” However the derivation in [4] of the QC bracket intentionally
avoids any kind of semiclassic limits due to its potential danger, see such an
attempt in [2] and Example 2 below. The actual method evaluates the image of
the universal bracket (2) under the QC representation [4, (20)] of the group
$\mathbb{D}^{m}$.

The paper [1, 3001-p4] “corrected” the original derivation of QC bracket
replacing the initial set of Planck constants $h_{1}$ and to $h_{2}$ by the
new one $h_{\mathrm{eff}}$ defined by the expression:

$\frac{1}{h_{\mathrm{eff}}}=\frac{1}{h_{1}}+\frac{1}{h_{2}}.$ (7)

However this transformation is singular for $h_{1}h_{2}=0$ and needs special
clarifications how to proceed for such values.

######  Example 2.

Let us consider the transformation $U_{h}:f(x,y)\mapsto f(hx,\frac{1}{h}y)$,
which is a unitary operator $L_{2}(\mathbb{R}^{2})\rightarrow
L_{2}(\mathbb{R}^{2})$ for any $h>0$. However this does not allow us “to take
the limit $h\rightarrow 0$” through the straightforward substitution $h=0$.

Furthermore the paper [1, 3001-p3] claims that “we have shown that the
equation of motion (16) does not lead to a non-trivial QC limit”. However,
this is caused by p-mechanisation (3)–(4), cf. the next two examples.

######  Example 3.

Let $B_{1}$ and $B_{2}$ are squares of coordinate $Q$ and momentum $P$
observables (of the quantum subsystem) respectively. Under p-mechanisation
(3)–(4) they are represented by squares of corresponding convolutions. Then
the commutator (first term of bracket (1)) of their QC representations is
zero, the second term in (1)) vanishes since no classical observables present,
and the third termin (1)) is equal to QC image of the observable $4QP$. Thus
the total bracket is indeed the same as the Poisson bracket for those
observables.

Let us examine the above claim for the p-mechanisation (5) and assume that two
p-mechanical observables $B_{1}$ and $B_{2}$, that is two convolutions on the
group $\mathbb{D}^{n}{}$ [4, p. 876], for any fixed $g_{1}$ are multiples of
the delta function in $g_{2}$, e.g. as in Ex. 3. Under the QC representation
$\rho_{(h;q,p)}$ [4, (20)] those observables become operators
$\rho_{(h;q,p)}(B_{1})$ and $\rho_{(h;q,p)}(B_{2})$ on the state space for the
quantum subsystem without any dependence from classical coordinates $p$, $q$
and the respective Planck constant $h_{2}$. Correspondingly the second and the
third terms of the bracket (1) vanish and this bracket is equal to the
(quantum) commutator
$\frac{1}{\mathrm{i}h}[\rho_{(h;q,p)}(B_{1}),\rho_{(h;q,p)}(B_{2})]$.

Therefore if we admit the claim [1, 3001-p3] that QC bracket (1) always
coincides with the purely classic Poisson bracket, then we have to accept that
any quantum commutator is always equal to the Poisson bracket.

######  Example 4.

Under p-mechanisation (5) the squares of momentum and coordinates from Ex. 3
are represented by convolutions with kernels
$\delta^{\prime\prime}_{x_{1}x_{1}}(g_{1};g_{2})$ and
$\delta^{\prime\prime}_{y_{1}y_{1}}(g_{1};g_{2})$. Their commutator on the
group $\mathbb{D}^{1}{}$ is
$4\delta^{\prime\prime\prime}_{x_{1}y_{1}s_{1}}+2\delta^{\prime\prime}_{s_{1}s_{1}}$.
Thus the universal bracket (2) is

$\left\\{\\!\left[B_{1},B_{2}\right]\\!\right\\}=4\delta^{\prime\prime}_{x_{1}y_{1}}+2\delta^{\prime}_{s_{1}}+(4\delta^{\prime\prime\prime}_{x_{1}y_{1}s_{1}}+2\delta^{\prime\prime}_{s_{1}s_{1}})\mathcal{A}_{2}.$
(8)

In the QC representation of $\mathbb{D}^{1}{}$ the last term of the sum
vanishes and two first terms produce $4QP+2\mathrm{i}hI$. This is the quantum
commutator of $Q^{2}$ and $P^{2}$ times $\frac{1}{\mathrm{i}h}$. There is no
unitary representation to get rid of the purely imaginary term $2\mathrm{i}hI$
in order to reduce the QC bracket of $B_{1}$ and $B_{2}$ to the value $4QP$ of
their Poisson bracket.

##  4 Conclusion

In this paper we demonstrated that the QC bracket (1) does not possess itself
two properties of “artificial coupling” and “genuinely classical nature” as
claimed in [1]. Unfortunately the claims [1] were uncritically translated by
some other authors, see [5, 6].

We showed that for a decoupled Hamiltonian the dynamics of observables
localised in one subsystem is unaffected by the Hamiltonian of the other
subsystem. The “classical nature” described in [1] is rooted in the
p-mechanisation used in [4] and does not appear with other choice of
p-mechanical observables.

The main conclusion of the commented paper [1] is: “We suggest that a
different Ansatz for the equations of motion, could indeed produce non-trivial
QC equations”. This comment is aimed to clarify possible directions for such a
search.

##  5 Acknowledgements

I am grateful to Dr. Frederica Agostini for useful discussions and providing
me with a part of her unpublished thesis. Prof. O.V. Prezhdo and an anonymous
referee made suggestions, which improved presentation in this comment.

## References

  * [1] Agostini F., Caprara S. Ciccotti G. Europhys. Lett. EPL 782007Art. 30001, 6 10.1209/0295-5075/78/30001.
  * [2] Prezhdo O. V. Kisil V. V. Phys. Rev. A (3) 561997162 arXiv:quant-ph/9610016.
  * [3] Kisil V. V. J. Phys. A 372004183 arXiv:quant-ph/0212101, On-line. .
  * [4] Kisil V. V. Europhys. Lett. 722005873 arXiv:quant-ph/0506122, On-line.
  * [5] Hall M. J. W. Physical Review A782008042104.
  * [6] Zhan F., Lin Y. Wu B. Journal of Chemical Physics 1282008315204\.


# Stellar spots cause measurable variations in atmospheric metallicity

Tanner A. Wilson,1,2 and Andrew R. Casey,1,2
1School of Physics & Astronomy, Monash University, Victoria, Australia
2Center of Excellence for Astrophysics in Three Dimensions (ASTRO-3D)
Corresponding Author E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

To accurately measure a star’s atmospheric parameters and chemical abundances,
it is crucial to have high-quality spectra. Analysing the detailed chemical
abundances of groups of stars can help us better understand nucleosynthesis,
galactic chemical enrichment, and stellar evolution. In this study, we
explored whether stellar spots can affect a star’s inferred metallicity and,
if so, where the impact is the strongest. To investigate this, we created
synthetic infrared spectra that included stellar spots for a sample of main-
sequence stars younger than the sun. We then applied two models to the data:
one that accounted for spots and one that did not. From this, we can determine
the bias introduced when fitting spotted spectra with a non-spotted model and
how this bias varies with different parameters. Our findings revealed that
fitting spotted spectra with a non-spotted model can introduce a scatter of up
to 0.05 dex in the inferred metallicity, especially for stars with high levels
of spot coverage. This bias is similar in magnitude to other relevant effects,
such as atomic diffusion, radiative levitation, or non-local thermodynamic
equilibrium. We also found that the effect is most pronounced in young stars
and decreases with age. These results suggest that stellar spots can introduce
a systematic uncertainty in metallicity that is not currently accounted for in
spectroscopic analysis. This could potentially limit scientific inferences for
population-level studies of young stars and differential abundance analyses.

###### keywords:

stars: abundances – stars:starspots – stars: rotation

††pubyear: 2023††pagerange: Stellar spots cause measurable variations in
atmospheric metallicity–References

##  1 Introduction

It is widely assumed that the elemental abundances in a star’s atmosphere
accurately reflect the abundances of the material from which the star formed
(Gibson et al., 2003; Pagel, 2009; Salaris & Cassisi, 2017). This assumption
is critical for chemical tagging (Anders et al., 2016; Randich et al., 2022),
understanding galactic formation (Gibson et al., 2003), and the synthesis of
elements across cosmic time (McWilliam & Rauch, 2004; Johnson et al., 2020).
Precise measurements of elemental abundances are essential in many areas of
astrophysics. For example, chemical tagging allows us to track the history of
the galaxy, which would be impossible with biased measures of abundance.
Differential abundance techniques (Önehag et al., 2011; Meléndez et al., 2014;
Reggiani et al., 2016; Maia et al., 2019; Liu et al., 2020; Nissen et al.,
2020; Spina et al., 2021) – employed for solar twins and planet-hosting stars
– claim very precise abundance measurements, which are essential for probing
planet formation (Tayar et al., 2022). Similarly, when determining cluster
ages (Bensby et al., 2004; Pont & Eyer, 2004), the turn-off age of a star is
particularly useful for this purpose because a small change in
colour/magnitude, which depends on metallicity, indicates a relatively large
change in age compared to the main-sequence.

Recognising that surface abundances may change over a star’s evolution is
important. The surface abundances can change due to numerous processes. Atomic
diffusion and radiative levitation introduce surface abundance variations on
the scale of 0.05 dex, with a magnitude and bias that depends on the element
and the stellar effective temperature (Önehag et al., 2014). Enhanced mixing
can also cycle material to the surface. Nuclear reactions, such as lithium
depletion (Pinsonneault et al., 2002) or CNO cycling (Crowther, 2007)) enhance
and deplete specific surface abundances and isotopic ratios. Accretion can
enhance surface metallicity and vary particular elemental abundances for a
short time depending on the companion type (Pasquini et al., 2007; Maldonado
et al., 2019; Laughlin & Adams, 1997). For example, mass loss can strip away
H-rich surface regions in Wolf-Rayet stars (Crowther, 2007) \- increasing the
observed stellar metallicity or carrying away surface metals which will have a
small to negligible decrease of surface metals on the main sequence.

These effects are usually ignored when estimating a star’s stellar parameters
and chemical abundances. Most spectroscopic analyses usually adopt some
simplifying assumptions to make the computation time tractable. For example,
we usually assume the stellar photosphere can be represented in one dimension
(1D) and that baryonic matter can be described by thermal distributions in
small regions (local thermal equilibrium; LTE). These assumptions can
particularly influence the measured stellar parameters (e.g., Blanco-Cuaresma,
2019). Both can lead to an over-estimate of the temperature gradient in the
atmosphere and an underestimation of the density, which can result in an over-
estimate of the abundance. We also typically ignore magnetic activity, but
recently Spina et al. (2020) showed it has a measurable impact on the chemical
abundances of young, fast-rotating stars. While these assumptions simplify
inference, it is important to consider their effects when reaching
conclusions.

Stars have spots, which are important indicators of the rotational rate of
stars, especially along the main sequence (McQuillan et al., 2014; Santos et
al., 2021). The properties of stellar spots and their effect on the observed
properties of a star vary with: age, rotation rate, mass, and metallicity
(Mathur et al., 2014; Karoff et al., 2018; Nichols-Fleming & Blackman, 2020).
For example, as the rotation rates of stars decrease with age, the average
magnetic activity likewise decreases. This results in smaller short lived
spots that cover only a small fraction of the stellar surface (Cao &
Pinsonneault, 2022). Spot properties can be generalised by: their coverage
across the stellar surface, the temperature difference relative to the
surrounding, and the occurrence pattern. Cao & Pinsonneault (2022) recently
quantified the spot parameters of stars in the Pleiades and M67. They found
that young or fast rotating stars tend to be more magnetically active and have
a greater spot coverage than their older, slower counterparts.

The spotted areas of the star can be thousands of degrees cooler than the
surrounding areas - solar spots for example can be 500-2000K cooler than the
surrounding photosphere (Berdyugina, 2005; Herbst et al., 2020). A spotted
star’s stellar spectra are more complex than their non-spotted counterparts
(Morris et al., 2019). Accurate inference of stellar parameters requires a
model that reflects the stellar spectra well. In this work, we quantify the
effect of fitting spotted spectra with non-spotted models and identify the
parts of the main sequence where the effect is most prevalent. In Section 2,
we outline the generative model for stellar spectra with spots and describe
our choices of stellar parameters before outlining the fitting procedure used.
In Section 3, we present the difference in the recovered stellar parameters
with the spotted and non-spotted models. We discuss parts of the main sequence
where the effect is most prevalent. Finally, in Section 4, we place those
results in the context of other significant effects on measured stellar
metallicity and provide recommendations for high-precision spectroscopic
investigations in specific regions of stellar evolution.

##  2 Method

###  2.1 Stellar parameters for a population of fake stars

We prepare a sample of stellar spectra that spans the main sequence to
estimate the impact that stellar spots can have on the accuracy of inferred
stellar parameters. This sample is intended to be indicative of a possible
population of main-sequence stars but not intended to represent which stars
would, or would not, have spots. We generate 1500 spectra of main-sequence and
early post-main-sequence stars with various values of mass, age, metallicity,
$v\sin{i}$, $f_{\text{spot}}$, and $x_{\text{spot}}$ across the HR diagram. We
drew masses from a Salpeter initial mass function (Salpeter, 1955) between 0.5
and 1.5 $M_{\odot}$ with $\alpha=2.35$. This limits our range of masses to
those with a radiative surface and convective core and reaches beyond the
Kraft break (Kraft, 1967). Metallicity is drawn from a distribution to
approximately reflect what is observed in the Milky Way. Specifically, we
defined a variable $\phi$ to be drawn from a Beta distribution

$\phi\sim\mathcal{B}\left(\alpha=10,\beta=2\right)$ (1)

and applied a transform from $\phi$ to [Fe/H] by requiring the metallicities
be bounded between $[\mathrm{Fe/H}]_{\mathrm{min}}=-2$ and
$[\mathrm{Fe/H}]_{\mathrm{max}}=+0.5$. We also required that the mode of
$\phi$, defined as $\frac{\alpha-1}{\alpha+\beta-2}$ for a Beta distribution,
occurs at Solar metallicity. This leads to the transform:

$[\mathrm{Fe/H}]=\left(\frac{}{}[\mathrm{Fe/H}]_{\mathrm{max}}-[\mathrm{Fe/H}]_{\mathrm{min}}\right)\left(\phi-\frac{\alpha-1}{\alpha+\beta-2}\right)\quad.$
(2)

The stars we generate mock data for in this work span from the zero-age main
sequence (ZAMS) to low-luminosity subgiants. We draw equivalent evolutionary
phase (EEP) values from a uniform distribution EEP $\sim\mathcal{U}(200,450)$,
where $\mathcal{U}\left(x,y\right)$ denotes a uniform prior between x and y.
The bounds of this range (200 and 450) represent the ZAMS and the low-
luminosity subgiant phase, respectively. Using the EEP, mass and metallicity,
we interpolate a position along the MIST stellar isochrones (Morton, 2015) to
calculate the expected $T_{\text{eff}}$ and $\log{g}$ for each random star. We
also obtain the star’s age (post-ZAMS) that we can use in conjunction with the
other stellar parameters to determine rotational properties (see below). We
have limited the age of the stars we consider in this work up to the age of
the Sun. This is the range available for rotational rate and convective
turnover timescales from the sources we draw from in this work. This limits
the post-MS stars we consider to more massive stars. We briefly discuss bias’
which may introduce in Section 4.

The surface rotation period is interpolated from stellar cluster-tuned
rotational isochrones given the stellar age and mass (Table A1 in Spada et al.
(2016)). Rotational broadening $v\sin{i}$ can then be calculated by combining
the rotational period, the radius from the interpolated isochrone model, and
an inclination angle. We have drawn inclination from a uniform distribution in
$\cos{i}\sim\mathcal{U}(0,1)$.

$f_{\text{spot}}$ is related to the Rossby number, $R_{o}$, which is defined
as the ratio of the surface rotational period to the convective turnover
timescale ($\tau_{\textrm{conv}}$). $\tau_{\textrm{conv}}$ is interpolated
from Table 1 in Landin et al. (2010) given the stellar age and mass. Combining
this value with the rotational period, we obtain $R_{o}$. $f_{\text{spot}}$ is
then calculated from the relationship between $f_{\text{spot}}$ and $Ro$
identified in Cao & Pinsonneault (2022) (Eq. 5):

$f_{\text{spot}}=\left\\{\begin{array}[]{lr}0.278,&\log{R_{o}}\leq-0.729\\\
0.0623\ R_{o}^{-0.881},&\log{R_{o}}>-0.729\end{array}\right\\}.$ (3)

There is some scatter in $f_{\text{spot}}$ which is not accounted for by this
relation (see left panel of Figure 7 in Cao & Pinsonneault (2022)). For this
reason, we add random noise to our calculated $f_{\text{spot}}$ which is drawn
from a normal distribution with a standard deviation of 0.1 it is clear that
there is some scatter in $f_{\text{spot}}$ About

We assume $x_{\text{spot}}$ is drawn from a uniform distribution
$x_{\text{spot}}\sim\mathcal{U}(0.8,1.0)$. This represents the limits set when
fitting $x_{\text{spot}}$ in (Cao & Pinsonneault, 2022), which is motivated by
temperature bounds which they discuss in more detail in Section 2.2.
$x_{\text{spot}}$ does not appear to have a clear relationship with other
stellar parameters, but it - and $f_{\text{spot}}$ \- may vary on multiple
periodic timescales as they do for the Sun. The stochastic nature of stellar
observations - and the admittedly simple nature of the model - means that
$f_{\text{spot}}$ and $x_{\text{spot}}$ are random draws from the possible
stellar spot parameters. We will eventually find that $x_{\text{spot}}$ has
little effect on the bias introduced by fitting spotted spectra with a non-
spotted model, so move forward with the knowledge that we have good coverage
when modelling over the range of possible parameters.

Figure 1: HR diagram of the 1500 sets of stellar parameters drawn from
physically motivated distributions of mass, metallicity and age coloured by
[Fe/H].

###  2.2 Spotted spectrum generative model

We build upon the work of (Cao & Pinsonneault, 2022), where a forward model is
developed to model the effect of starspots and to estimate the fractional spot
coverage of stars in the Pleiades and M67. Their model assumes that the
spectrum of a spotted star can be broken into spotted and non-spotted
components. These two components have the same $\log{g}$ [Fe/H],
microturbulent velocity, and the same surface rotational velocity
($v\sin{i}$), but the two components vary in temperature. The spot and ambient
temperatures ($T_{\text{spot}}$ and $T_{\text{amb}}$) are related by
$T_{\text{spot}}$ = $x_{\text{spot}}$$T_{\text{amb}}$, and are coupled to the
effective temperature of the star following the approach of Somers &
Pinsonneault (2015)

$T_{\text{eff}}=T_{\text{amb}}(1-f_{\text{spot}}+f_{\text{spot}}x_{\text{spot}}^{4})^{\frac{1}{4}},$
(4)

where $f_{\text{spot}}$is the fractional surface spot coverage. From these
relations, the set {$T_{\text{eff}}$, $x_{\text{spot}}$,$f_{\text{spot}}$}
define a pair of ambient and spot temperatures that preserve stellar
luminosity.

We calculated a grid of synthetic spectra, which we interpolate between to
generate the predicted spectra for a spotted or non-spotted model. The list of
atomic and molecular transitions is from (Shetrone et al., 2015; Smith et al.,
2021). We used a grid of plane-parallel MARCS (Gustafsson et al., 2008) model
photospheres that span dimensions in effective temperature, surface gravity,
and metallicity.111We calculated spectra using spherical models as well, but
in practice, only spectra from plane-parallel models (i.e., main-sequence
stars) are used in this work. Microturbulence was kept fixed at 1.15 km s-1
for main-sequence stars and we assumed that [$\alpha$/H] scales with [Fe/H]
(i.e., the so-called ‘standard’ composition in MARCS). The abundance
dimensions [C/M] and [N/M] were kept fixed at zero. We used Korg (Wheeler et
al., 2022) to synthesise all model spectra at high resolution, which we then
convolved and down-sampled to match the (uniform in log) pixel spacing used in
the APOGEE data reduction pipeline (Holtzman et al., 2018). The convolution
kernel includes two components that enter multiplicatively: one assuming a
constant spectral resolution $R=\lambda/\Delta\lambda$ of 22,500, and another
representing rotational broadening $v\sin{i}$. We convolved each spectrum with
a grid of $v\sin{i}$ values that were uniformly spaced in $\log{v\sin{i}}$
from 0-100 km s-1 in order to match the setup for the APOGEE analysis
pipeline. Naturally, for low $v\sin{i}$ values, the line spread function of
the instrument will dominate.

With this grid of spectra and some given spectral parameters
$\\{T_{\mathrm{eff}},\log{g},[\mathrm{Fe/H}],\log{v\sin{i}}\\}$, we
interpolate the spotted and ambient spectra and combine them in a fractional
manner with wavelength as if they were separate black-body spectra in order to
produce a flux-preserving combined spectrum:

$B(T_{\text{eff}},\lambda)\ =\
f_{\text{spot}}B(T_{\text{spot}},\lambda)+(1-f_{\text{spot}})B(T_{\text{amb}},\lambda)\quad.$
(5)

In total, our forward model for predicting spotted spectra includes six
parameters: $T_{\mathrm{eff}}$, $\log{g}$, [Fe/H], $\log{v\sin{i}}$,
$x_{\text{spot}}$, and $f_{\text{spot}}$. This model is equally capable of
predicting non-spotted spectra by fixing $f_{\text{spot}}$ to zero or
$x_{\text{spot}}$ to unity.

Using the 1500 sets of parameters outlined in Section 2.1 we generated
synthetic spotted stellar spectra. We also apply realistic noise at each pixel
from a Gaussian distribution with standard deviation = 0.01, assuming a
signal-to-noise ratio of 100. Continuum normalisation is performed by assuming
a running mean of the spectra, and during fitting, this procedure is applied
to the fake spectrum (data) and to the model spectrum.

We now have the tools to determine the effect of fitting spotted spectra with
non-spotted models. We do this by finding the best-fitting stellar parameters
given the synthetic spectra fitted twice: first with the model described in
Section 2.2 and then with a non-spotted model (e.g., $f_{\text{spot}}$ fixed
at 0 and $x_{\text{spot}}$ fixed at 1). Here we have performed least-squares
fitting through the Levenberg-Marquardt algorithm implemented in SciPy. We
found that fitting the spotted parameters can be non-trivial. The likelihood
surfaces are multimodal and degenerate, requiring informed choices about the
initialisation of fitting. To resolve this issue, we performed a coarse
evaluation of parameters (on a grid) before starting optimisation.

##  3 Results

We began by confirming that we could accurately recover the injected
parameters. The best-fit parameters following fitting the synthetic spotted
spectra with a spotted model are shown in Figures 2 and 3. $T_{\text{eff}}$,
$\log{g}$, [Fe/H], and $v\sin{i}$ are recovered accurately for every injected
parameter set. While we identify scatter in recovered $f_{\text{spot}}$ this
appears not to affect the accuracy of the recovery of the traditional stellar
parameters. We move forward confident that any difference in the recovered
parameters between fitting with the spotted and non-spotted models results
from the model differences rather than the fitting procedure employed in this
work.

Figure 2: Recovered traditional stellar parameters ($T_{\text{eff}}$,
$\log{g}$, [Fe/H] and $v\sin{i}$) from fitting synthetic spotted spectra with
a spotted model of the stellar atmosphere against the corresponding injected
parameters. We consistently accurately recover each injected value when a
spotted model of the stellar atmosphere is employed to fit the spotted
synthetic spectra. Figure 3: Recovered spot parameters ($x_{\text{spot}}$ and
$f_{\text{spot}}$) from the synthetic spotted spectra fitted with a spotted
model of the stellar spectra against the injected parameters of the synthetic
spectra. We identify that the spot parameters are not always accurately
recovered through the fitting procedure. The recovered spot parameters are
notably more inaccurate as $x_{\text{spot}}$ approaches 1.  Figure 4: The
difference between the recovered traditional stellar spectra parameters
($T_{\text{eff}}$, $\log{g}$, [Fe/H] and $v\sin{i}$) from the synthetic
spotted spectra fitted with both a spotted and non-spotted model of the
stellar spectra against the injected parameters of the synthetic spectra
(spotted model non-spotted model recovered parameter). We identify scatter
introduced to each of the stellar parameters when fitting spotted spectra with
a non-spotted model of the stellar atmosphere.

We now identify systematic effects in the recovered parameters when we fit the
spotted spectra with an incorrect model of non-spotted spectra. The difference
between the recovered parameters fitted with a spotted and non-spotted model
of the stellar atmosphere are shown in Figure 4. A consistent scatter is
introduced on each parameter when a non-spotted model is used to perform
inference on a spotted spectrum. We calculate each parameter’s average bias
and scatter to quantify the effect. The injected parameters are separated into
10 bins, and we take the median and median absolute deviation of the
difference between the spotted and non-spotted model’s inferred parameters for
each bin. We take the median as a measure of the average bias and the median
absolute deviation as a proxy for the scatter.

In Figure 5 we show the effect of the injected parameters on the stellar spot
spectra through the difference between the recovered spot and non-spot model
$T_{\text{eff}}$. Fitting a spotted spectrum with a small
$x_{\text{spot}}$with a non-spotted spectrum introduces a consistent bias to
the inferred $T_{\text{eff}}$ of about $-25$ K: a non-spotted model tends to
underestimate the true effective temperature of a spotted spectrum. A scatter
is also introduced $T_{\text{eff}}$ on the scale of $\sim$50K for spectra with
significant spot coverage (low $x_{\text{spot}}$ and large $f_{\text{spot}}$).
The other injected parameters do not appear to have any strong correlations or
effects on the recovered non-spot $T_{\text{eff}}$. Their median values are
zero, and MAD appears consistent at $\sim$25K.

Figure 5: Bias introduced to $T_{\text{eff}}$ (blue) when fitting spotted
spectra with a non-spotted model against injected parameters of synthetic
spectra. Each injected parameter is binned into ten bins across the range of
injected parameters. The median and median absolute deviation of the
difference between the spotted and non-spotted recovered $T_{\text{eff}}$
($\Delta$$T_{\text{eff}}$) are then calculated in each bin. Scatter points
show the median $\Delta$$T_{\text{eff}}$ for each bin in injected parameters.
Filled areas show one maximum absolute deviation above and below the median
value and dashed lines indicate the edge of this range. Inference of
$T_{\text{eff}}$ with a non-spotted model injects random scatter on average of
the scale of $\sim$50K and introduces a consistent bias of order $\sim$-25K
for spectra with significant spot coverage.  Figure 6: Bias introduced to
$\log{g}$ (orange), [Fe/H] (green) and $\log$ $v\sin{i}$ (red) when fitting
spotted spectra with a non-spotted model against injected parameters of
synthetic spectra. Here we have plotted $\log$ $v\sin{i}$ as the grid we
perform interpolation over is distributed uniformly in $\log$ $v\sin{i}$ and
to make direct comparisons to the other stellar parameters ([Fe/H] and
$\log{g}$). Each injected parameter is binned into ten bins across the range
of injected parameters. The median and median absolute deviation of the
difference between the spotted and non-spotted recovered parameters ($\Delta$
parameter) are then calculated in each of these bins. Scatter points show each
bin’s median $\Delta$ parameter in injected parameters. Filled areas show one
maximum absolute deviation above and below the median value and dashed lines
indicate the edge of this range. Fitting spotted spectra with non-spotted
models does not introduce a bias to the inferred parameters though it does
introduce a scatter. This scatter increases with decreasing $x_{\text{spot}}$
and increasing $f_{\text{spot}}$ to a maximum median absolute deviation on the
order of $\sim$0.05 for $\log{g}$, [Fe/H] and $\log$ $v\sin{i}$ \- which
corresponds to a maximum scatter on $v\sin{i}$ of $\sim$2 km/s. The scatter is
otherwise approximately constant for all other parameters, and on average
$\sim 0.025$ \- which corresponds to an average scatter on $v\sin{i}$ of
$\sim$ 1 km/s

Figure 6 shows the effects of fitting spotted spectra with a non-spotted model
on $\log{g}$ (orange) and $\log$ $v\sin{i}$ (red), respectively. There appears
to be no statistically significant bias introduced to both of the inferred
parameters as the median of each bin of injected parameters is consistently
about zero. However, a consistent scatter is introduced to both parameters.
The MAD of $\Delta$ $\log{g}$ and $\Delta\log$ $v\sin{i}$ in each injected
parameter bin have an average value of $\sim$0.025 dex - corresponding to an
average scatter on $v\sin{i}$ of $\sim$ 1 km s-1. The scatter peaks for both
recovered parameters at $\sim$0.05 dex for stars with significant spot
coverage - which corresponds to a maximum scatter on $v\sin{i}$ of $\sim$ 2
km/s. The scatter on recovered $v\sin{i}$ and $\log{g}$ is otherwise constant
with the other injected parameters.

The effect of fitting spotted spectra with a non-spotted model is significant
in the recovery of metallicity. This is seen in Figure 6 (green), where we
compare the recovered [Fe/H] with a spotted and non-spotted model of the
stellar atmosphere against the injected parameters of our spotted spectra.
This process does not introduce a bias to the inferred metallicity of the
spectra but does introduce a significant scatter to the recovered value,
representing an intrinsic ‘minimum floor’ of systematic uncertainty if the
effects of spots are not included (see Section 4).

The scatter introduced to [Fe/H] by fitting spotted spectra with a non-spotted
model increases with injected $f_{\text{spot}}$. As $f_{\text{spot}}$
approaches 1, the MAD of $\Delta$[Fe/H] reaches a maximum of about 0.05 dex.
Comparatively, as $x_{\text{spot}}$ decreases, so does the MAD of
$\Delta$[Fe/H], peaking again at 0.04 dex. As the scatter in the other
injected parameters is relatively constant, there is no significant relation
between the other spectral parameters and $\Delta$[Fe/H]. The introduced
scatter in [Fe/H] is dominated by the spot parameters of spectra.

Figure 7: Bias introduced to [Fe/H] when fitting spotted spectra with a non-
spotted model against the age of the model used to generate synthetic spectra.
Age is binned into ten bins across the range of injected parameters. The
median and median absolute deviation of the difference between the spotted and
non-spotted recovered [Fe/H] ($\Delta$[Fe/H]) are then calculated in each of
these bins. Scatter points show the median $\Delta$[Fe/H] for each bin in
injected parameters. Filled areas show one maximum absolute deviation above
and below the median value and dashed lines indicate the edge of this range.
We find that the introduced scatter is greatest for stars younger than $\sim$
2Gyr is $\sim$0.02 while the bias is 0 for stars in this age range. The
scatter decreases for older stars ($>$3Gyr) to $\sim$0.01 but the median
$\Delta$[Fe/H] increases with increasing age. This is like the result of the
small number of stars in our sample in this age range rather than indicative
of a trend.

##  4 Discussion

The results in Section 3 indicate that using a non-spotted model to fit
spotted stellar spectra introduces a systematic bias of up to $-25$ K in
effective temperature and no substantial bias in other parameters. In their
study of fitting a spotted stellar model to APOGEE spectra of members in the
Pleiades and M67, Cao & Pinsonneault (2022) find a systematic 0.1 dex
enhancement in observed [Fe/H]. The lack of bias we find here could be
attributed to different stellar populations of stars (e.g., some stars are
biased in one direction, but in our population, that effect is mitigated by
biases in the opposite direction). We find that the effects of model mismatch
(i.e., using a non-spotted model to fit a spotted spectrum) can also introduce
a scatter (measured by median absolute deviation) of about 50 K in effective
temperature and 0.05 dex in other parameters. If we assume that the spot model
we adopt is representative of reality, then these scatter values would
represent a minimum systematic uncertainty in these parameters if the wrong
model (a non-spotted model) is used. These deviations are comparable to the
typical random uncertainties reported by the APOGEE survey (150K, 0.13 dex and
0.1 dex; Hegedűs et al., 2022), although these random uncertainties will vary
with signal-to-noise.

Systematic uncertainties (like model mismatches) will dominate in high signal-
to-noise ratios, and the level of scatter we find in metallicity (0.05 dex) is
comparable to the effects of radiative levitation, atomic diffusion (Önehag et
al., 2014), and magnetic broadening of absorption lines (Spina et al., 2020).
Unlike these effects, which can in part be mitigated through parameterisation
with other stellar parameters, accounting for stellar spots requires a model
that explicitly predicts their contribution to the emergent spectrum. This
scatter in [Fe/H] is significant as it is of the same order as the precision
of spectroscopic inference of metallicity. In particular, a differential
analysis of two Solar twins might report abundance uncertainties on the level
of 0.01-0.02 dex. While the two stars are selected to be extremely similar in
order to mitigate systematic effects, those two stars could have very
different coverages of stellar spots, which would introduce a systematic
uncertainty floor.

###  4.1 Imperfect models

The results we show here are limited in their applicability. When generating
the mock data, only a fraction of randomly drawn stellar parameters could be
used to synthesise spectra, either because of limitations of stellar
isochrones, the spectral grid, or limits in the procedure in estimating an
appropriate rotational velocity and Rossby number. We also limit the ages of
the stellar sample to 4.6 Gyr - the maximum ages of both the models used to
determine the convective turnover timescale and grid of rotational periods set
by observations. As a result, our sample is limited to relatively young stars,
and there are hints of a bias in injected parameters towards higher
$f_{\text{spot}}$. We have extensively probed the region where the effect
should be most prevalent in terms of the scatter it introduces, but this is
not intended to be a complete and representative population of main-sequence
stars. The quantitative results may not be perfectly accurate for some regions
of the HR diagram, and might vary with photopshere geometry. However, by
assuming spots are present everywhere across the main sequence, our analysis
shows where the consequential effects are most or least prevalent.

The treatment of stellar spots in this work requires some discussion. Stellar
spots are highly complex regions on the surface of stars. The position of
spots relative to the observer, their temporal evolution, and the inherent
magnetic activity and faculae surrounding stellar spots, would all introduce
complexity to the emergent spectra from these regions. The spotted model
employed in this work is a first-order approximation of the average effect of
spots on stellar spectra. The functional form of the temporal evolution of the
stellar spots in stars other than the Sun is not well known. For a given
$f_{\text{spot}}$ we could assume that $x_{\text{spot}}$ varies on some
periodic or temporal scale, even if we don’t know the functional form of that
variability. In this scenario with our model, $x_{\text{spot}}$ is drawn from
a uniform prior, which implicitly assumes that we are observing the star at
some random time. This modelling of $x_{\text{spot}}$ is relatively crude
since, in principle, $x_{\text{spot}}$ could vary as a function of other
stellar parameters.

Investigations of the evolution of fractional spot coverage of stars is a
developing field. For example, recent works have shown an enhancement in
$f_{\text{spot}}$ for stars undergoing core-envelope recoupling (Cao et al.,
2023). For this reason, our results are only indicative rather than
prescriptive. Applying this model to more stars APOGEE samples and time series
spectroscopic observations of stars could elucidate the relationship between
the parameters.

Cao & Pinsonneault (2022) suggest that young, magnetically active stars -
stars with Rossby numbers < 0.4 - have $f_{\text{spot}}$ greater than > 0.1,
saturating at $f_{\text{spot}}$ $\sim$ 0.3, with significant scatter, when
$R_{o}$ < 0.2. There is also a significant scatter in $f_{\text{spot}}$ for
these stars. The use of the Rossby number to reflect the magnetic/spot
activity of stars should be treated with some care. The Sun expresses periodic
evolution of its magnetic activity (time scale on the order of decades) and
stellar spot expression (time scale on the order of years). The range of
fractional spot coverages we observe in the Sun is on the order of [0, 0.12]
without variations in the Rossby number. As a result, we draw the injected
$f_{\text{spot}}$from relations with Rossby number and add a random scatter
drawn from a Gaussian distribution with a standard deviation of 0.1.

Employing a non-spotted spectra model to fit spotted spectra can introduce a
scatter to recovered parameters, but fitting a spotted model to non-spotted
spectra has little to no effect on the recovered parameters. We recommend that
a spotted model, if only as simple as the one used in this work, will
consistently recover stellar parameters better than a non-spotted model while
also providing a measure of the spot parameters of stars.

###  4.2 When should a spotted model of the stellar atmosphere be employed?

The scatter introduced to the recovered stellar parameters increases with
fractional spot coverage. Fractional spot coverage is inversely related to the
rotation rate of stars through $R_{o}$. Further, the rotation rate of stars
decreases with time, owing to magnetic braking. As a result, the fractional
spot coverage of stars is expected to decrease with age.

We can probe when the scatter introduced to the recovery of stellar parameters
by stellar spots is most prominent by calculating the scatter in
$\Delta$[Fe/H] with bins of age. In Figure 7, we show the bias and scatter
introduced to $\Delta$[Fe/H] with respect to stellar age. The introduced
scatter is greatest for stars younger than $\sim$ 2Gyr is $\sim$0.02, while
the bias, measured through the median, is zero for stars in this age range.
The scatter decreases for older stars ($>$3Gyr) to $\sim$0.01, but the median
$\Delta$[Fe/H] appears to increase with increasing age. The increase in the
median value is most likely not indicative of a trend and rather the result of
the low number of stars in the larger age bins.

The trends that we identify in this work are only qualitative - though they do
allow us to make recommendations for future work. Our results indicate that
fitting the spotted spectra of a star with a non-spotted model when
$f_{\text{spot}}$> 0.1 will introduce a scatter to bias the recovered
parameters. We suggest using a spotted model if a star is significantly
photometrically variable due to stellar spots. McQuillan et al. (2014)
calculated the rotation periods of low-mass main-sequence stars that are
photometrically variable due to stellar spots. They were able to determine the
rotation rates of stars across a wide mass range (0.6 < $M_{\odot}$ < 1.1) at
multiple points along the main sequence. These stars must therefore express
stellar spots and may have the measured stellar parameters influenced by the
effect we identify in this work. They estimated that $\sim$23% of main-
sequence stars exhibit definite rotational modulation from stellar spots, a
lower bound due to observational effects. We, therefore, believe at least 1/4
of the main sequence stars may be affected by this bias.

##  5 Conclusions

Here we have shown that stellar spots can introduce measurable systematic bias
and variance to inferred stellar parameters when a non-spotted model is used.
The results demonstrate that spectra with strong spot features can introduce a
scatter in inferred metallicity of order 0.05 dex. This emphasises the need
for caution when performing spectroscopic analysis on stars with visible
spots, particularly young, fast-rotating stars. Our findings highlight the
importance of incorporating the effect of spots into spectroscopic models to
ensure accurate and precise results.

The magnitude of this effect is comparable to others that plague stellar
spectroscopy, including atomic diffusion, radiative levitation, and non-local
thermodynamic equilibrium. However, the impact of this effect will vary
depending on the scientific context. Turn-off ages of clusters are likely to
be only minimally impacted, as the metallicity bias for old, slowly rotating
stars is less than 0.01 dex. In contrast, a systematic error floor of 0.05 dex
caused by spots on the main sequence would critically limit the capacity of
strong chemical tagging (Casamiquela et al., 2021). Similarly, star spots
could limit any inferences from differential abundance analyses of Sun-like
stars, where the typical reported uncertainty is 0.01-0.02 dex (e.g., Meléndez
et al., 2014; Nissen, 2015; Reggiani et al., 2016; Maia et al., 2019; Liu et
al., 2020; Nissen et al., 2020; Spina et al., 2021). While we have focused on
the impact on overall metallicity and not on individual abundances, it will be
important to examine these effects more closely at a per-element level. These
results provide valuable insights for future studies on stars and their
properties and underscore the need for continued research on the impact of
spots on spectroscopic inference.

## Acknowledgements

We thank the reviewer for their constructive feedback which helped improve the
paper. We thank Prof. Ilya Mandel for his helpful feedback and our fruitful
discussions. A. R. C. thanks Jon Holtzman and Adam Wheeler for helpful
discussions. A. R. C. is supported in part by the Australian Research Council
through a Discovery Early Career Researcher Award (DE190100656). Parts of this
research were supported by the Australian Research Council Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project
number CE170100013.

## Data Availability

The data and models underlying this article are available upon request to the
corresponding author.

## References

  * Anders et al. (2016) Anders F., et al., 2016, Astronomische Nachrichten, 337, 926
  * Bensby et al. (2004) Bensby T., Feltzing S., Lundström I., 2004, Astronomy & Astrophysics, 421, 969
  * Berdyugina (2005) Berdyugina S. V., 2005, Living Reviews in Solar Physics, 2, 8
  * Blanco-Cuaresma (2019) Blanco-Cuaresma S., 2019, Monthly Notices of the Royal Astronomical Society, 486, 2075
  * Cao & Pinsonneault (2022) Cao L., Pinsonneault M. H., 2022, Monthly Notices of the Royal Astronomical Society, 517, 2165
  * Cao et al. (2023) Cao L., Pinsonneault M. H., van Saders J. L., 2023, Core-envelope decoupling drives radial shear dynamos in cool stars, http://arxiv.org/abs/2301.07716
  * Casamiquela et al. (2021) Casamiquela L., Castro-Ginard A., Anders F., Soubiran C., 2021, Astronomy & Astrophysics, 654, A151
  * Crowther (2007) Crowther P. A., 2007, Annual Review of Astronomy & Astrophysics, vol. 45, Issue 1, pp.177-219, 45, 177
  * Gibson et al. (2003) Gibson B. K., Fenner Y., Renda A., Kawata D., Lee H., 2003, Publications of the Astronomical Society of Australia, 20, 393
  * Gustafsson et al. (2008) Gustafsson B., Edvardsson B., Eriksson K., Jørgensen U. G., Nordlund \\., Plez B., 2008, Astronomy & Astrophysics, 486, 951
  * Hegedűs et al. (2022) Hegedűs V., Mészáros S., Jofré P., Stringfellow G. S., Feuillet D., García-Hernández D. A., Nitschelm C., Zamora O., 2022, Comparative analysis of atmospheric parameters from high-resolution spectroscopic sky surveys: APOGEE, GALAH, Gaia-ESO, doi:10.48550/arXiv.2211.03416, https://arxiv.org/abs/2211.03416v2
  * Herbst et al. (2020) Herbst K., Papaioannou A., Airapetian V., Atri D., 2020, The Astrophysical Journal, 907
  * Holtzman et al. (2018) Holtzman J. A., et al., 2018, The Astronomical Journal, 156, 125
  * Johnson et al. (2020) Johnson J. A., Fields B. D., Thompson T. A., 2020, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 378, 20190301
  * Karoff et al. (2018) Karoff C., et al., 2018, The Astrophysical Journal, 852, 46
  * Kraft (1967) Kraft R. P., 1967, The Astrophysical Journal, 150, 551
  * Landin et al. (2010) Landin N. R., Mendes L. T. S., Vaz L. P. R., 2010, Astronomy and Astrophysics, 510, A46
  * Laughlin & Adams (1997) Laughlin G., Adams F. C., 1997, The Astrophysical Journal, 491, L51
  * Liu et al. (2020) Liu F., Yong D., Asplund M., Wang H. S., Spina L., Acuña L., Meléndez J., Ramírez I., 2020, Monthly Notices of the Royal Astronomical Society, 495, 3961
  * Maia et al. (2019) Maia M. T., Meléndez J., Lorenzo-Oliveira D., Spina L., Jofré P., 2019, Astronomy and Astrophysics, 628, A126
  * Maldonado et al. (2019) Maldonado J., Villaver E., Eiroa C., Micela G., 2019, Astronomy & Astrophysics, 624, A94
  * Mathur et al. (2014) Mathur S., et al., 2014, Astronomy and Astrophysics, 562, A124
  * McQuillan et al. (2014) McQuillan A., Mazeh T., Aigrain S., 2014, The Astrophysical Journal Supplement Series, 211, 24
  * McWilliam & Rauch (2004) McWilliam A., Rauch M., eds, 2004, Origin and evolution of the elements. No. v. 4 in Carnegie Observatories astrophysics series, Cambridge University Press, Cambridge ; New York
  * Meléndez et al. (2014) Meléndez J., et al., 2014, The Astrophysical Journal, 791, 14
  * Morris et al. (2019) Morris B. M., Curtis J. L., Sakari C., Hawley S. L., Agol E., 2019, The Astronomical Journal, 158, 101
  * Morton (2015) Morton T. D., 2015, Astrophysics Source Code Library, p. ascl:1503.010
  * Nichols-Fleming & Blackman (2020) Nichols-Fleming F., Blackman E., 2020, Monthly Notices of the Royal Astronomical Society, 491, 2706
  * Nissen (2015) Nissen P. E., 2015, Astronomy & Astrophysics, 579, A52
  * Nissen et al. (2020) Nissen P. E., Christensen-Dalsgaard J., Mosumgaard J. R., Silva Aguirre V., Spitoni E., Verma K., 2020, Astronomy and Astrophysics, 640, A81
  * Pagel (2009) Pagel B., 2009, Nucleosynthesis and Chemical Evolution of Galaxies, by Bernard E. J. Pagel, Cambridge, UK: Cambridge University Press, 2009
  * Pasquini et al. (2007) Pasquini L., Döllinger M. P., Weiss A., Girardi L., Chavero C., Hatzes A. P., Silva L. d., Setiawan J., 2007, Astronomy & Astrophysics, 473, 979
  * Pinsonneault et al. (2002) Pinsonneault M. H., Steigman G., Walker T. P., Narayanan V. K., 2002, The Astrophysical Journal, 574, 398
  * Pont & Eyer (2004) Pont F., Eyer L., 2004, Monthly Notices of the Royal Astronomical Society, 351, 487
  * Randich et al. (2022) Randich S., et al., 2022, A&A, 666
  * Reggiani et al. (2016) Reggiani H., Meléndez J., Yong D., Ramírez I., Asplund M., 2016, Astronomy & Astrophysics, 586, A67
  * Salaris & Cassisi (2017) Salaris M., Cassisi S., 2017, Royal Society Open Science, 4, 170192
  * Salpeter (1955) Salpeter E. E., 1955, The Astrophysical Journal, 121, 161
  * Santos et al. (2021) Santos A. R. G., Breton S. N., Mathur S., García R. A., 2021, The Astrophysical Journal Supplement Series, 255, 17
  * Shetrone et al. (2015) Shetrone M., et al., 2015, The Astrophysical Journal Supplement Series, 221, 24
  * Smith et al. (2021) Smith V. V., et al., 2021, The Astronomical Journal, 161, 254
  * Somers & Pinsonneault (2015) Somers G., Pinsonneault M. H., 2015, Proceedings of the International Astronomical Union, 10, 91
  * Spada et al. (2016) Spada F., Gellert M., Arlt R., Deheuvels S., 2016, Astronomy & Astrophysics, 589, A23
  * Spina et al. (2020) Spina L., et al., 2020, The Astrophysical Journal, 895, 52
  * Spina et al. (2021) Spina L., Sharma P., Meléndez J., Bedell M., Casey A. R., Carlos M., Franciosini E., Vallenari A., 2021, Nature Astronomy, 5, 1163
  * Tayar et al. (2022) Tayar J., Claytor Z. R., Huber D., van Saders J., 2022, The Astrophysical Journal, 927, 31
  * Wheeler et al. (2022) Wheeler A. J., Abruzzo M. W., Casey A. R., Ness M. K., 2022, Korg: a modern 1D LTE spectral synthesis package, http://arxiv.org/abs/2211.00029
  * Önehag et al. (2011) Önehag A., Korn A., Gustafsson B., Stempels E., Vandenberg D. A., 2011, Astronomy and Astrophysics, 528, A85
  * Önehag et al. (2014) Önehag A., Gustafsson B., Korn A., 2014, Astronomy & Astrophysics, 562, A102