00footnotetext:  The authors are supported by the grant
FRGS/1/2019/STG06/UM/02/6

# Geodesic deviation equation in $f(Q)$-gravity

Jing-Theng Beh and Tee-How Loo and Avik De  J. T. Beh
Institute of Mathematical Sciences
Universiti Malaya
50603 Kuala Lumpur
Malaysia<EMAIL_ADDRESS>T. H. Loo
Institute of Mathematical Sciences
Universiti Malaya
50603 Kuala Lumpur
Malaysia<EMAIL_ADDRESS>A. De
Department of Mathematical and Actuarial Sciences
Universiti Tunku Abdul Rahman
Jalan Sungai Long
43000 Cheras
Malaysia<EMAIL_ADDRESS>
###### Abstract.

In the present paper we study the Geodesic Deviation Equation (GDE) in the
modified $f(Q)$-gravity theories. The formulation of GDE in General Relativity
in the case of the homogeneous and isotropic Friedman-Lemaître-Robertson-
Walker (FLRW) spacetime is briefly discussed and then extended in modified
$f(Q)$-gravity using its covariant counterpart. The generalised Mattig
relation is obtained. Finally, an equivalent expression to the Dyer-Roeder
equation in General Relativity in the case of $f(Q)$-gravity is presented.

##  1\. INTRODUCTION

Einstein’s General Relativity (GR) is one of the most successful theories in
Physics. This set of simple looking but complicated enough field equations

$R_{\mu\nu}-\frac{R}{2}g_{\mu\nu}=\kappa T_{\mu\nu}.$ (1)

has provided a remarkable narrative of the cosmological observational data
[19], and created new insights into the concepts of space and time. The
mathematical framework of this geometrical theory of gravity is based on
(pseudo-)Riemannian geometry. It describes the properties of the gravitational
field by using the curvature tensor of the spacetime. One of the fundamental
equations in this theory is the geodesic deviation equation (GDE) which
provides a relation between the Riemannian curvature tensor and the relative
acceleration between two nearby test particles. This equation describes the
relative motion of free falling particles to bend towards or away from each
other, under the impact of a gravitational field [18]. However, despite this
undeniable success, increasing technological ability in modern observational
cosmology posed new questions to GR. It turns out that the validity of GR
might only be up to the astrophysical scales not exceeding the Solar system
[14, 3].

To resolve the imperfection of GR, one of the approaches is to modify the
matter sector of the field equations (1) by adding some additional ‘dark’
components to the energy budget of the universe, and the other one is to
modify the gravitational sector. The most common modifications in the latter
direction are achieved by generalizing the Einstein-Hilbert action term,
precisely by replacing the Ricci scalar $R$ with an arbitrary function of some
scalar produced from the Riemann or Ricci curvature tensor or some topological
invariant, like the Gauss-Bonnet term or combination of both scalar-tensor
terms; $f(R)$ gravity theory being the simplest and most popular one such [16,
7]. However, GR is formulated based on a very special and unique connection
coefficient, the torsionless and metric-compatible Levi-Civita connection
which can be posed as a function of the metric tensor and apparently there are
other gravity theories equivalent to GR, such as the teleparallel gravity (TG)
[10] and symmetric teleparallel gravity (STG) [13] unhindered to this special
connection coefficient. Unlike GR, where the gravity is described by the
curvature of spacetime, in both these theories the curvature is set to be
zero, and the torsion of the connection controls the gravity in TG and the
non-metricity of the connection does so in STG. As theories equivalent to GR,
these two also inherit the same ‘dark’ problem and so following the same idea
to consider the dark contents of the universe as the contribution of the
spacetime geometry as in $f(R)$-theory, modification of these two theories was
also due. In this way, the role of Ricci scalar $R$ in GR is replaced by a
scalar $T$ formed from the torsion tensor $T^{\alpha}_{\>\>\beta\gamma}$ in TG
[11] and by scalar $Q$ formed from the non-metricity tensor
$Q_{\alpha\beta\gamma}$ in STG [12] and thus modified $f(T)$ and $f(Q)$
theories were born. Both these theories have some drawback coined as “the good
and bad tetrad problem” in TG [17] and “the good and bad coordinates problem”
in STG [21]. Like in TG, where the consistency of the theory depends on the
choice of tetrad, the $f(Q)$-theory, which is the main focus of the present
article, also relies on the choice of coordinates. In fact, under the
constraint of vanishing curvature and torsion, we can always choose the so-
called ‘coincident gauge’ in which the affine connection vanishes and take
metric as the only fundamental variable, however, the theory will no longer be
diffeomorphism invariant and might be inconsistent in some coordinates system.
To avoid this issue, $f(Q)$ theory can be formulated in a covariant way [21].

As a natural extension to GR, the GDE was formulated in $f(R)$-gravity [9, 8].
Although the notion of geodesic deviation in TG is slightly different than GR
in the sense that in GR, the motion of particle is described by the curvature
of spacetime, so the GDE serves as the force equation, and on the other hand,
in TG, the torsion appeared as a real force, namely, the tidal force; the
teleparallel depiction of the gravitational interaction is totally equivalent
to that of GR [1]. Thus, it is completely natural to put the force equation in
TG into the form of geodesic equation in GR. In this approach, the
corresponding GDE in TG can be obtained [5]. This motivated us to investigate
the covariant formulation of STG and find the GDE in the $f(Q)$-theory.

The outline of this paper is as follows: After introduction, in section 2 we
reformulate the covariant version of field equations of $f(Q)$-theory. In
section 3, we recapitulate the GDE in GR briefly before we plunge into the
same in $f(Q)$-theory in section 4. Next we consider the ansatz of RW metric
in section 5 and then further discuss the case of fundamental observer and
null vector fields in this setting in section 6 and section 7, respectively.
We finish with a discussion of Dyer-Roeder like equation in section 8 and a
concluding section at the end.

##  2\. FIELD EQUATIONS

We begin with constructing a non-metric affine connection in symmetric
teleparallelism, that is, $\nabla_{\lambda}g_{\mu\nu}\neq 0$ and we define the
non-metricity tensor

$Q_{\lambda\mu\nu}:=\nabla_{\lambda}g_{\mu\nu}\,,$ (2)

so $Q_{\lambda\mu\nu}=Q_{\lambda(\mu\nu)}$. The associated affine connection
can be expressed as

$\Gamma^{\lambda}{}_{\mu\nu}:=\mathring{\Gamma}^{\lambda}{}_{\mu\nu}+L^{\lambda}{}_{\mu\nu}$
(3)

where $\mathring{\Gamma}^{\lambda}{}_{\mu\nu}$ is the Levi-Civita connection
from the metric

$\mathring{\Gamma}^{\lambda}{}_{\mu\nu}=\frac{1}{2}g^{\lambda\rho}(\partial_{\mu}g_{\rho\nu}+\partial_{\nu}g_{\mu\rho}-\partial_{\rho}g_{\mu\nu})$

and $L^{\lambda}{}_{\mu\nu}$ is called the disformation tensor. It follows
that,

$L^{\lambda}{}_{\mu\nu}=\frac{1}{2}(Q^{\lambda}{}_{\mu\nu}-Q_{\mu}{}^{\lambda}{}_{\nu}-Q_{\nu}{}^{\lambda}{}_{\mu})\,.$
(4)

From the definition of non-metricity tensor, we construct two different types
of non-metricity vectors,

$Q_{\mu}:=g^{\nu\lambda}Q_{\mu\nu\lambda}=Q_{\mu}{}^{\nu}{}_{\nu}\,,\qquad\tilde{Q}_{\mu}:=g^{\nu\lambda}Q_{\nu\mu\lambda}=Q_{\nu\mu}{}^{\nu}\,.$

Then, from (4), we have

$L^{\alpha}{}_{\alpha\mu}=-\frac{1}{2}Q_{\mu}\,,\qquad
L_{\mu\alpha}{}^{\alpha}=\frac{1}{2}Q_{\mu}-\tilde{Q}_{\mu}\,.$

Moreover, we define the trace of the non-metricity tensor

$Q:=g^{\mu\nu}(L^{\alpha}{}_{\beta\mu}L^{\beta}{}_{\nu\alpha}-L^{\alpha}{}_{\beta\alpha}L^{\beta}{}_{\mu\nu})$
(5)

and the superpotential tensor

$P^{\lambda}{}_{\mu\nu}:=\frac{1}{4}\left(-2L^{\lambda}{}_{\mu\nu}+Q^{\lambda}g_{\mu\nu}-\tilde{Q}^{\lambda}g_{\mu\nu}-\frac{1}{2}\delta^{\lambda}_{\mu}Q_{\nu}-\frac{1}{2}\delta^{\lambda}_{\nu}Q_{\mu}\right)\,.$
(6)

Therefore, from (4), (5) and (6), we have [20]

$Q=Q_{\lambda\mu\nu}P^{\lambda\mu\nu}=-\frac{1}{2}Q_{\lambda\mu\nu}L^{\lambda\mu\nu}+\frac{1}{4}Q_{\lambda}Q^{\lambda}-\frac{1}{2}Q_{\lambda}\tilde{Q}^{\lambda}\,.$
(7)

As discussed in [20], the $f(Q)$-theory was constructed by using the
constraints, $R^{\rho}{}_{\sigma\mu\nu}=0$. That means there exist a special
coordinate system such that the affine connection vanishes,
$\Gamma^{\lambda}{}_{\mu\nu}=0$. This situation is called the coincident
gauge. Under this circumstance, the metric is the only dynamical variable. But
as mentioned in [21], in any other coordinate such that the connection does
not vanish, the evolution of the metric will be affected, and results a
completely different theory. Therefore, by varying the action term

$S=\int\left[\frac{1}{2\kappa}f(Q)+\mathcal{L}_{M}\right]\sqrt{-g}\,d^{4}x$

with respect to thmetric, where $\kappa=8\pi G$, we obtain the field equations
[20]

$\frac{2}{\sqrt{-g}}\nabla_{\lambda}(\sqrt{-g}f^{\prime}P^{\lambda}{}_{\mu\nu})-\frac{1}{2}fg_{\mu\nu}+f^{\prime}(P_{\nu\rho\sigma}Q_{\mu}{}^{\rho\sigma}-2P_{\rho\sigma\mu}Q^{\rho\sigma}{}_{\nu})=\kappa
T_{\mu\nu}$ (8)

where $f^{\prime}=\partial f/\partial Q$. However, this equation is only valid
in the coincident gauge coordinate.

On the other hand, we know that the curvature tensor can be written as

$R^{\rho}{}_{\sigma\mu\nu}=\partial_{\mu}\Gamma^{\rho}{}_{\nu\sigma}-\partial_{\nu}\Gamma^{\rho}{}_{\mu\sigma}+\Gamma^{\rho}{}_{\mu\lambda}\Gamma^{\lambda}{}_{\nu\sigma}-\Gamma^{\rho}{}_{\nu\lambda}\Gamma^{\lambda}{}_{\mu\sigma}\,.$

Then, by using (3), we found that

$R^{\rho}{}_{\sigma\mu\nu}=\mathring{R}^{\rho}{}_{\sigma\mu\nu}+\mathring{\nabla}_{\mu}L^{\rho}{}_{\nu\sigma}-\mathring{\nabla}_{\nu}L^{\rho}{}_{\mu\sigma}+L^{\rho}{}_{\mu\lambda}L^{\lambda}{}_{\nu\sigma}-L^{\rho}{}_{\nu\lambda}L^{\lambda}{}_{\mu\sigma}$

and so

$\displaystyle R_{\sigma\nu}$
$\displaystyle=\mathring{R}_{\sigma\nu}+\frac{1}{2}\mathring{\nabla}_{\nu}Q_{\sigma}+\mathring{\nabla}_{\rho}L^{\rho}{}_{\nu\sigma}-\frac{1}{2}Q_{\lambda}L^{\lambda}{}_{\nu\sigma}-L^{\rho}{}_{\nu\lambda}L^{\lambda}{}_{\rho\sigma}$
$\displaystyle R$
$\displaystyle=\mathring{R}+\mathring{\nabla}_{\lambda}Q^{\lambda}-\mathring{\nabla}_{\lambda}\tilde{Q}^{\lambda}-\frac{1}{4}Q_{\lambda}Q^{\lambda}+\frac{1}{2}Q_{\lambda}\tilde{Q}^{\lambda}-L_{\rho\nu\lambda}L^{\lambda\rho\nu}\,.$

Thus, from the teleparallelism, $R^{\rho}{}_{\sigma\mu\nu}=0$, we obtain

$\mathring{R}_{\mu\nu}-\frac{1}{2}\mathring{R}=2\nabla_{\lambda}P^{\lambda}_{\mu\nu}-\frac{1}{2}Qg_{\mu\nu}+(P_{\rho\mu\nu}Q^{\rho\sigma}{}_{\sigma}+P_{\nu\rho\sigma}Q_{\mu}{}^{\rho\sigma}-2P_{\rho\sigma\mu}Q^{\rho\sigma}{}_{\nu})\,.$

It follows that the field equations in (8) can be rewritten as [21]

$f^{\prime}\mathring{G}_{\mu\nu}+\frac{1}{2}g_{\mu\nu}(f^{\prime}Q-f)+2f^{\prime\prime}\nabla_{\lambda}QP^{\lambda}{}_{\mu\nu}=\kappa
T_{\mu\nu}$ (9)

where

$\mathring{G}_{\mu\nu}=\mathring{R}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\mathring{R}$

and $T_{\mu\nu}$ is the energy-momentum tensor. This equation is in a form
similar to the field equations in $f(R)$-gravity and we can proceed to find
the GDE.

##  3\. GEODESIC DEVIATION EQUATION IN GR

In this section, we depict several concepts about the geodesic deviation
equation (GDE) in GR. Consider a congruence of geodesics of arbitrary casual
character described by $x^{\alpha}(\nu,s)$, where $\nu$ is an affine parameter
of geodesics, and $s$ is an indexed family of geodesics. That is, for each
fixed $s$, the curve described by $\gamma_{s}(\nu)$ is a geodesic. Let
$V^{\alpha}$ denote the normalized tangent vector field of the congruence,
then $V^{\alpha}=\frac{dx^{\alpha}}{d\nu}$ and
$V_{\alpha}V^{\alpha}=\epsilon$, where $\epsilon=+1,0,-1$ if the geodesics are
spacelike, null, or timelike respectively. Define
$\eta^{\alpha}=\frac{dx^{\alpha}}{ds}$ as the deviation vector for the
congruence.

Since $V^{\alpha}$ and $\eta^{\alpha}$ commutes, that is,
$\mathcal{L}_{V}\eta^{\alpha}=\mathcal{L}_{\eta}V^{\alpha}$ (or
$[V,\eta]^{\alpha}=0$), so
$\nabla_{V}\nabla_{V}\eta^{\alpha}=\nabla_{V}\nabla_{\eta}V^{\alpha}$. Then,
by using the Ricci identity,
$\nabla_{X}\nabla_{Y}Z^{\alpha}-\nabla_{Y}\nabla_{X}Z^{\alpha}-\nabla_{[X,Y]}Z^{\alpha}=\mathring{R}^{\alpha}{}_{\beta\gamma\delta}Z^{\beta}X^{\gamma}Y^{\delta}$,
we obtain the GDE [18]

$\frac{D^{2}\eta^{\alpha}}{D\nu^{2}}=-\mathring{R}^{\alpha}{}_{\beta\gamma\delta}V^{\beta}\eta^{\gamma}V^{\delta}$
(10)

where $\frac{D}{D\nu}$ is the covariant derivative along the geodesic.

To further simplify, we assume the energy-momentum tensor in the form of a
perfect fluid

$T_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu}$ (11)

where $\rho$ is the energy density and $p$ is the pressure. It follows that
the trace of the energy-momentum tensor is

$T=3p-\rho\,.$ (12)

Recall that we have the Einstein field equations in GR (with cosmological
constant)

$\mathring{R}_{\mu\nu}-\frac{1}{2}\mathring{R}g_{\mu\nu}+\Lambda
g_{\mu\nu}=\kappa T_{\mu\nu}\,.$

Then, by using (11) and (12), the Ricci scalar and Ricci tensor can be
expressed as

$\mathring{R}=\kappa(\rho-3p)+4\Lambda$
$\mathring{R}_{\mu\nu}=\kappa(\rho+p)u_{\mu}u_{\nu}+\frac{1}{2}[\kappa(\rho-p)+2\Lambda]g_{\mu\nu}\,.$

Moreover, the Riemannian curvature tensor can also be expressed as [18]

$\displaystyle\mathring{R}_{\alpha\beta\gamma\delta}=C_{\alpha\beta\gamma\delta}$
$\displaystyle+\frac{1}{2}(g_{\alpha\gamma}\mathring{R}_{\delta\beta}-g_{\alpha\delta}\mathring{R}_{\gamma\beta}+g_{\beta\delta}\mathring{R}_{\gamma\alpha}-g_{\beta\gamma}\mathring{R}_{\delta\alpha})$
$\displaystyle-\frac{\mathring{R}}{6}\left(g_{\alpha\gamma}g_{\delta\beta}-g_{\alpha\delta}g_{\gamma\beta}\right)$
(13)

where $C_{\alpha\beta\gamma\delta}$ is the Weyl tensor. If we consider
$C_{\alpha\beta\gamma\delta}=0$, together with
$\epsilon=V_{\alpha}V^{\alpha},E=-V_{\alpha}u^{\alpha}$, and
$\eta_{\alpha}V^{\alpha}=\eta_{\alpha}u^{\alpha}=0$, then the right hand side
of GDE in (10) can be simplified as

$\mathring{R}^{\alpha}{}_{\beta\gamma\delta}V^{\beta}\eta^{\gamma}V^{\delta}=\bigg{[}\frac{1}{3}(\kappa\rho+\Lambda)\epsilon+\frac{1}{2}\kappa(\rho+p)E^{2}\bigg{]}\eta^{\alpha}\,.$

This is the well-known Pirani equation [6].

##  4\. GEODESIC DEVIATION EQUATION IN $f(Q)$-GRAVITY

In this section, we formulate the GDE in $f(Q)$-gravity. By contracting the
field equations in (9) with $g_{\mu\nu}$, we obtain the Ricci scalar

$\mathring{R}=\frac{1}{f^{\prime}}(2f^{\prime}Q-2f+2f^{\prime\prime}P^{\lambda\rho}{}_{\rho}\nabla_{\lambda}Q-\kappa
T)\,.$

Then, substituting this Ricci scalar into (9), we have the Ricci tensor

$\displaystyle\mathring{R}_{\mu\nu}=\frac{1}{f^{\prime}}\bigg{[}$
$\displaystyle\frac{1}{2}g_{\mu\nu}(f^{\prime}Q-f+2f^{\prime\prime}P^{\lambda\rho}{}_{\rho}\nabla_{\lambda}Q-\kappa
T)$
$\displaystyle-2f^{\prime\prime}P^{\lambda}{}_{\mu\nu}\nabla_{\lambda}Q+\kappa
T_{\mu\nu}\bigg{]}\,.$

Hence, by using (3) and considering $C_{\alpha\beta\gamma\delta}=0$, we found
that

$\displaystyle\mathring{R}_{\alpha\beta\gamma\delta}=$
$\displaystyle\frac{1}{2f^{\prime}}\bigg{[}\kappa(g_{\alpha\gamma}T_{\delta\beta}-g_{\alpha\delta}T_{\gamma\beta}+g_{\beta\delta}T_{\gamma\alpha}-g_{\beta\gamma}T_{\delta\alpha})$
$\displaystyle+\left(\frac{f^{\prime}Q}{3}-\frac{f}{3}-\frac{2\kappa
T}{3}+\frac{4}{3}f^{\prime\prime}P^{\lambda\rho}{}_{\rho}\nabla_{\lambda}Q\right)(g_{\alpha\gamma}g_{\delta\beta}-g_{\alpha\delta}g_{\gamma\beta})$
$\displaystyle+(g_{\alpha\gamma}\mathcal{D}_{\delta\beta}-g_{\alpha\delta}\mathcal{D}_{\gamma\beta}+g_{\beta\delta}\mathcal{D}_{\gamma\alpha}-g_{\beta\gamma}\mathcal{D}_{\delta\alpha})f^{\prime}\bigg{]}$

where

$\mathcal{D}_{\mu\nu}:=-2P^{\lambda}{}_{\mu\nu}\nabla_{\lambda}Q\,\partial_{Q}\,.$
(14)

Assume the perfect fluid form of the energy-momentum tensor stated in (11) and
(12), the above equation reduces to

$\displaystyle\mathring{R}_{\alpha\beta\gamma\delta}=$
$\displaystyle\frac{1}{2f^{\prime}}\bigg{[}\kappa(\rho+p)(g_{\alpha\gamma}u_{\delta}u_{\beta}-g_{\alpha\delta}u_{\gamma}u_{\beta}+g_{\beta\delta}u_{\gamma}u_{\alpha}-g_{\beta\gamma}u_{\delta}u_{\alpha})$
$\displaystyle+\left(\frac{f^{\prime}Q}{3}-\frac{f}{3}+\frac{2\kappa\rho}{3}+\frac{4}{3}f^{\prime\prime}P^{\lambda\rho}{}_{\rho}\nabla_{\lambda}Q\right)(g_{\alpha\gamma}g_{\delta\beta}-g_{\alpha\delta}g_{\gamma\beta})$
$\displaystyle+(g_{\alpha\gamma}\mathcal{D}_{\delta\beta}-g_{\alpha\delta}\mathcal{D}_{\gamma\beta}+g_{\beta\delta}\mathcal{D}_{\gamma\alpha}-g_{\beta\gamma}\mathcal{D}_{\delta\alpha})f^{\prime}\bigg{]}\,.$

Then, contracting with $V^{\beta}V^{\delta}$ and consider
$V_{\alpha}V^{\alpha}=\epsilon$, we have

$\displaystyle\mathring{R}_{\alpha\beta\gamma\delta}V^{\beta}V^{\delta}=$
$\displaystyle\frac{1}{2f^{\prime}}\bigg{[}\kappa(\rho+p)[g_{\alpha\gamma}(u_{\beta}V^{\beta})^{2}-2(u_{\beta}V^{\beta})V_{(\alpha}u_{\gamma)}+\epsilon
u_{\alpha}u_{\gamma}]$
$\displaystyle+\left(\frac{f^{\prime}Q}{3}-\frac{f}{3}+\frac{2\kappa\rho}{3}+\frac{4}{3}f^{\prime\prime}P^{\lambda\rho}{}_{\rho}\nabla_{\lambda}Q\right)(\epsilon
g_{\alpha\gamma}-V_{\alpha}V_{\gamma})$
$\displaystyle+[(g_{\alpha\gamma}\mathcal{D}_{\delta\beta}-g_{\alpha\delta}\mathcal{D}_{\gamma\beta}+g_{\beta\delta}\mathcal{D}_{\gamma\alpha}-g_{\beta\gamma}\mathcal{D}_{\delta\alpha})f^{\prime}]V^{\beta}V^{\delta}\bigg{]}\,.$

By raising the $\alpha$ index in the Riemannian tensor and contracting with
$\eta^{\gamma}$, we obtain

$\displaystyle\mathring{R}^{\alpha}{}_{\beta\gamma\delta}V^{\beta}\eta^{\gamma}V^{\delta}=$
$\displaystyle\frac{1}{2f^{\prime}}\bigg{[}\kappa(\rho+p)[(u_{\beta}V^{\beta})^{2}\eta^{\alpha}-(u_{\beta}V^{\beta})V^{\alpha}(u_{\gamma}\eta^{\gamma})$
$\displaystyle-(u_{\beta}V^{\beta})u^{\alpha}(V_{\gamma}\eta^{\gamma})+\epsilon
u^{\alpha}(u_{\gamma}\eta^{\gamma})]$
$\displaystyle+\left(\frac{f^{\prime}Q}{3}-\frac{f}{3}+\frac{2\kappa\rho}{3}+\frac{4}{3}f^{\prime\prime}P^{\lambda\rho}{}_{\rho}\nabla_{\lambda}Q\right)(\epsilon\eta^{\alpha}-V^{\alpha}(V_{\gamma}\eta^{\gamma}))$
$\displaystyle+[(\delta_{\gamma}^{\alpha}\mathcal{D}_{\delta\beta}-\delta^{\alpha}_{\delta}\mathcal{D}_{\gamma\beta}+g_{\beta\delta}\mathcal{D}_{\gamma}^{\alpha}-g_{\beta\gamma}\mathcal{D}_{\delta}^{\alpha})f^{\prime}]V^{\beta}\eta^{\gamma}V^{\delta}\bigg{]}\,.$
(15)

By using $-V_{\alpha}u^{\alpha}=E$ and
$\eta_{\alpha}u^{\alpha}=\eta_{\alpha}V^{\alpha}=0$, (4) becomes

$\displaystyle\mathring{R}^{\alpha}{}_{\beta\gamma\delta}V^{\beta}\eta^{\gamma}V^{\delta}$
$\displaystyle=\frac{1}{2f^{\prime}}\bigg{[}\kappa(\rho+p)E^{2}+\epsilon\left(\frac{2\kappa\rho}{3}+\frac{f^{\prime}Q}{3}-\frac{f}{3}+\frac{4}{3}f^{\prime\prime}P^{\lambda\rho}{}_{\rho}\nabla_{\lambda}Q\right)\bigg{]}\eta^{\alpha}$
$\displaystyle+\frac{1}{2f^{\prime}}\bigg{[}(\delta^{\alpha}_{\gamma}\mathcal{D}_{\delta\beta}-\delta^{\alpha}_{\delta}\mathcal{D}_{\gamma\beta}+g_{\beta\delta}\mathcal{D}^{\alpha}_{\gamma}-g_{\beta\gamma}\mathcal{D}^{\alpha}_{\delta})f^{\prime}V^{\beta}V^{\delta}\bigg{]}\eta^{\gamma}\,.$
(16)

##  5\. GDE with FLRW background

Assuming that the universe is homogeneous and isotropic, and described by the
spatially flat Friedmann-Robertson-Walker (FLRW) metric, where the line
element in the Cartesian coordinates is given by

$ds^{2}=-dt^{2}+a^{2}(t)\delta_{ij}dx^{i}dx^{j}$ (17)

where $a(t)$ is the scale factor. This implies that the only non-vanishing
Christoffel symbols are [18]

$\mathring{\Gamma}^{l}{}_{0j}=\frac{\dot{a}}{a}\delta^{l}_{j}=\mathring{\Gamma}^{l}{}_{j0},\qquad\mathring{\Gamma}^{0}{}_{ij}=a\dot{a}\delta_{ij}$
(18)

here $i,j,k,...=1,2,3$. So the Ricci scalar can be written as

$\mathring{R}=6\frac{\ddot{a}}{a}+6\left(\frac{\dot{a}}{a}\right)^{2}\,.$ (19)

In addition, the FLRW metric is conformally flat, that is,
$C_{\alpha\beta\gamma\delta}=0$.

Alternatively, the Christoffel symbol can be expressed as

$\mathring{\Gamma}_{\lambda\mu\nu}=-\frac{\dot{a}}{a}(-u_{\lambda}g_{\mu\nu}+u_{\mu}g_{\nu\lambda}+u_{\nu}g_{\lambda\mu}+u_{\lambda}u_{\mu}u_{\nu})\,.$

It follows that by using (17), (18) and (19) can be easily obtained.

Moreover, by using the spatially flat FLRW metric in (17), we find that (from
appendix)

$Q=-6H^{2}$ (20)

where $H:=\frac{\dot{a}}{a}$ is the Hubble parameter. Therefore, we know that
$Q$ is only time-dependent, so

$\nabla_{\lambda}Q=12H\dot{H}u_{\lambda}\,.$ (21)

Then, by using (20), (21) and the operator defined in (14), we find that some
of the specific terms in (4) can be expressed as (from appendix)

$\displaystyle(\delta^{\alpha}_{\gamma}\mathcal{D}_{\delta\beta}-\delta^{\alpha}_{\delta}\mathcal{D}_{\gamma\beta}+g_{\beta\delta}\mathcal{D}^{\alpha}_{\gamma}-g_{\beta\gamma}\mathcal{D}^{\alpha}_{\delta})f^{\prime}V^{\beta}\eta^{\gamma}V^{\delta}=-24H^{2}\dot{H}f^{\prime\prime}(2\epsilon+E^{2})\eta^{\alpha}$
$\displaystyle\frac{4}{3}\epsilon\eta^{\alpha}f^{\prime\prime}P^{\lambda\rho}{}_{\rho}\nabla_{\lambda}Q=48H^{2}\dot{H}f^{\prime\prime}\epsilon\eta^{\alpha}\,.$

Therefore, (4) reduces to

$\displaystyle\mathring{R}^{\alpha}{}_{\beta\gamma\delta}V^{\beta}\eta^{\gamma}V^{\delta}=$
$\displaystyle\frac{1}{2f^{\prime}}\bigg{[}(\kappa\rho+\kappa
p-24H^{2}\dot{H}f^{\prime\prime})E^{2}$
$\displaystyle+\left(\frac{2\kappa\rho}{3}+\frac{f^{\prime}Q}{3}-\frac{f}{3}\right)\epsilon\bigg{]}\eta^{\alpha}$
(22)

which is considered to be the generalized Pirani equation. Finally, the GDE in
$f(Q)$ gravity can be written as

$\displaystyle\frac{D^{2}\eta^{\alpha}}{D\nu^{2}}=$
$\displaystyle-\frac{1}{2f^{\prime}}\bigg{[}(\kappa\rho+\kappa
p-24H^{2}\dot{H}f^{\prime\prime})E^{2}+\left(\frac{2\kappa\rho}{3}+\frac{f^{\prime}Q}{3}-\frac{f}{3}\right)\epsilon\bigg{]}\eta^{\alpha}\,.$

Notice that in this GDE only the magnitude of the deviation vector
$\eta^{\alpha}$ is changed along the geodesics, which reflects the homogeneity
and isotropy of the FLRW universe. This is not the case in anistropic
universes, such as Bianchi I, where the GDE also induces a change in the
direction of the deviation vector, as shown in [4].

##  6\. GDE for fundamental observers with FLRW background

In the case of fundamental observers, as the affine parameter, $\nu$ coincides
with the proper time of the fundamental observer, $t$, we have
$V^{\alpha}=u^{\alpha}$. This implies that $\epsilon=-1$ and $E=1$. Then, (5)
reduces to

$\mathring{R}^{\alpha}{}_{\beta\gamma\delta}u^{\beta}\eta^{\gamma}u^{\delta}=\frac{1}{f^{\prime}}\left(\frac{\kappa\rho}{6}+\frac{\kappa
p}{2}-\frac{f^{\prime}Q}{6}+\frac{f}{6}-12H^{2}\dot{H}f^{\prime\prime}\right)\eta^{\alpha}\,.$
(23)

If we let $\eta^{\alpha}=le^{\alpha}$, where $e^{\alpha}$ is parallel
propagated along $t$, then isotropy leads to

$\frac{De^{\alpha}}{Dt}=0$

and so

$\frac{D^{2}\eta^{\alpha}}{Dt^{2}}=\frac{d^{2}l}{dt^{2}}e^{\alpha}\,.$

Thus, by using (10) and (23), we obtain

$\frac{d^{2}l}{dt^{2}}=-\frac{1}{f^{\prime}}\left(\frac{\kappa\rho}{6}+\frac{\kappa
p}{2}-\frac{f^{\prime}Q}{6}+\frac{f}{6}-12H^{2}\dot{H}f^{\prime\prime}\right)l\,.$

By letting $l=a(t)$, we have

$\frac{\ddot{a}}{a}=-\frac{1}{f^{\prime}}\left(\frac{\kappa\rho}{6}+\frac{\kappa
p}{2}-\frac{f^{\prime}Q}{6}+\frac{f}{6}-12H^{2}\dot{H}f^{\prime\prime}\right)\,.$
(24)

This equation is a special case of the generalized Raychaudhuri equation.
Notice that the above equation can also be obtained by the standard forms of
the modified Friedmann equations in the $f(Q)$-gravity model for flat universe
[2]

$3H^{2}=\frac{1}{f^{\prime}}\left[\kappa\rho+\frac{1}{2}(f^{\prime}Q-f)\right]$
(25) $2\dot{H}+3H^{2}=-\frac{1}{f^{\prime}}\left[\kappa
p-\frac{1}{2}(f^{\prime}Q-f)-24H^{2}\dot{H}f^{\prime\prime}\right]\,.$

##  7\. GDE for null vector fields with FLRW background

In the case of past-directed null vector fields, we have
$V^{\alpha}=k^{\alpha}$ with $k_{\alpha}k^{\alpha}=0$, and so $\epsilon=0$.
Then, (5) becomes

$\mathring{R}^{\alpha}{}_{\beta\gamma\delta}k^{\beta}\eta^{\gamma}k^{\delta}=\frac{1}{2f^{\prime}}(\kappa\rho+\kappa
p-24H^{2}\dot{H}f^{\prime\prime})E^{2}\eta^{\alpha}\,.$ (26)

This equation can be explained as the Ricci focusing in $f(Q)$-gravity. If we
consider $\eta^{\alpha}=\eta
e^{\alpha},\,e_{\alpha}e^{\alpha}=1,\,e_{\alpha}u^{\alpha}=e_{\alpha}k^{\alpha}=0$
and $\frac{De^{\alpha}}{D\nu}=k^{\beta}\nabla_{\beta}e^{\alpha}=0$, in which
the basis is both aligned and propagated, then (5) can be written in a new
form

$\frac{d^{2}\eta}{d\nu^{2}}=-\frac{1}{2f^{\prime}}(\kappa\rho+\kappa
p-24H^{2}\dot{H}f^{\prime\prime})E^{2}\eta\,.$ (27)

As in the case of GR [6], all past-directed null geodesics experience focusing
if $\kappa(\rho+p)>0$ except the special case with the equation of state
$p=-\rho$. Thus, it is clear that (27) indicates the focusing condition for
the $f(Q)$-gravity, which is

$\frac{\kappa(\rho+p)}{f^{\prime}}>\frac{24H^{2}\dot{H}f^{\prime\prime}}{f^{\prime}}\,.$

After that, (27) can be expressed in term of redshift parameter, $z$. First,
we write

$\frac{d}{d\nu}=\frac{dz}{d\nu}\frac{d}{dz}$

which implies that

$\displaystyle\frac{d^{2}}{d\nu^{2}}$
$\displaystyle=\frac{dz}{d\nu}\frac{d}{dz}\left(\frac{d}{d\nu}\right)$
$\displaystyle=\left(\frac{d\nu}{dz}\right)^{-2}\left[-\left(\frac{d\nu}{dz}\right)^{-1}\frac{d^{2}\nu}{dz^{2}}\frac{d}{dz}+\frac{d^{2}}{dz^{2}}\right]\,.$
(28)

For the null geodesics, we have

$(1+z)=\frac{a_{0}}{a}=\frac{E}{E_{0}}\quad\longrightarrow\quad\frac{dz}{1+z}=-\frac{da}{a}$

where $a_{0}=1$ the present value of the scale factor. For the past-directed
case, we set $E_{0}=-1$, so

$dz=-(1+z)\frac{1}{a}\frac{da}{d\nu}d\nu=-(1+z)\frac{\dot{a}}{a}Ed\nu=H(1+z)^{2}d\nu$

and so

$\frac{d\nu}{dz}=\frac{1}{H(1+z)^{2}}$

Consequently,

$\frac{d^{2}\nu}{dz^{2}}=-\frac{1}{H(1+z)^{3}}\left[\frac{1}{H}(1+z)\frac{dH}{dz}+2\right]$

where

$\frac{dH}{dz}=\frac{d\nu}{dz}\frac{dt}{d\nu}\frac{dH}{dt}=-\frac{1}{H(1+z)}\frac{dH}{dt}$

where we make use of $\frac{dt}{d\nu}=E=-(1+z)$. Using
$\frac{\ddot{a}}{a}=\dot{H}+H^{2}$ in (24) we get

$\dot{H}=-\frac{1}{f^{\prime}}\left(\frac{\kappa\rho}{6}+\frac{\kappa
p}{2}-\frac{f^{\prime}Q}{6}+\frac{f}{6}-12H^{2}\dot{H}f^{\prime\prime}\right)-H^{2}\,.$

Hence,

$\displaystyle\frac{d^{2}\nu}{dz^{2}}=-\frac{3}{H(1+z)^{3}}\bigg{[}1+\frac{1}{3H^{2}f^{\prime}}\bigg{(}\frac{\kappa\rho}{6}+\frac{\kappa
p}{2}-$ $\displaystyle\frac{f^{\prime}Q}{6}+\frac{f}{6}$
$\displaystyle-12H^{2}\dot{H}f^{\prime\prime}\bigg{)}\bigg{]}\,.$

Putting this equation in (7), we have

$\displaystyle\frac{d^{2}\eta}{d\nu^{2}}=(H(1+z)^{2})^{2}\bigg{[}\frac{d^{2}\eta}{dz^{2}}+\frac{3}{(1+z)}$
$\displaystyle\bigg{[}1+\frac{1}{3H^{2}f^{\prime}}\bigg{(}\frac{\kappa\rho}{6}+\frac{\kappa
p}{2}$
$\displaystyle-\frac{f^{\prime}Q}{6}+\frac{f}{6}-12H^{2}\dot{H}f^{\prime\prime}\bigg{)}\bigg{]}\frac{d\eta}{dz}\bigg{]}.$

Finally, by using (27), the null GDE can be written in the form

$\displaystyle\frac{d^{2}\eta}{dz^{2}}$
$\displaystyle+\frac{3}{(1+z)}\left[1+\frac{1}{3H^{2}f^{\prime}}\left(\frac{\kappa\rho}{6}+\frac{\kappa
p}{2}-\frac{f^{\prime}Q}{6}+\frac{f}{6}-12H^{2}\dot{H}f^{\prime\prime}\right)\right]\frac{d\eta}{dz}$
$\displaystyle+\frac{\kappa(\rho+p)-24H^{2}\dot{H}f^{\prime\prime}}{2H^{2}(1+z)^{2}f^{\prime}}\eta=0\,.$
(29)

Since content of the universe is the ordinary matter and the radiation, so the
$\rho$ and $p$ can be be expressed as

$\kappa\rho=3H^{2}_{0}\Omega_{m_{0}}(1+z)^{3}+3H^{2}_{0}\Omega_{r_{0}}(1+z)^{4},\qquad\kappa
p=H^{2}_{0}\Omega_{r_{0}}(1+z)^{4}$ (30)

where we use $p_{m}=0$ and $p_{r}=\frac{1}{3}\rho_{r}$. From (25) and (30) ,
we could express $H^{2}$ as

$H^{2}=\frac{H^{2}_{0}}{f^{\prime}}[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{DE}]$
(31)

where

$\Omega_{DE}:=\frac{1}{H^{2}_{0}}\left(\frac{f^{\prime}Q}{6}-\frac{f}{6}\right)$
(32)

is the Dark Energy parameter. Hence, by using (30) and (31), the null GDE in
(7) can be expressed as

$\frac{d^{2}\eta}{dz^{2}}+\mathcal{P}(H,\dot{H},z)\frac{d\eta}{dz}+\mathcal{Q}(H,\dot{H},z)\eta=0$
(33)

where

$\displaystyle\mathcal{P}(H,\dot{H},z)=$
$\displaystyle\frac{\frac{7}{2}\Omega_{m_{0}}(1+z)^{3}+4\Omega_{r_{0}}(1+z)^{4}+2\Omega_{DE}}{(1+z)[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{DE}]}$
$\displaystyle-\frac{\frac{12\dot{H}f^{\prime\prime}}{f^{\prime}}[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{DE}]}{(1+z)[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{DE}]}$
(34) $\displaystyle\mathcal{Q}(H,\dot{H},z)=$
$\displaystyle\frac{3\Omega_{m_{0}}(1+z)^{3}+4\Omega_{r_{0}}(1+z)^{4}}{2(1+z)^{2}[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{DE}]}$
$\displaystyle-\frac{\frac{24\dot{H}f^{\prime\prime}}{f^{\prime}}[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{DE}]}{2(1+z)^{2}[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{DE}]}\,.$
(35)

In a particular case, where $f(Q)=Q-2\Lambda$, so $f^{\prime}=1$ and
$f^{\prime\prime}=0$. Thus, $\Omega_{DE}$ in (32) reduces to

$\Omega_{DE}=\frac{1}{H^{2}_{0}}\left(\frac{Q}{6}-\frac{Q-2\Lambda}{6}\right)=\frac{\Lambda}{3H^{2}_{0}}=:\Omega_{\Lambda}\,.$

This implies that the $H^{2}$ in (31) becomes the same as the Friedmann
equation in GR

$H^{2}=H^{2}_{0}[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{\Lambda}]$

which confirms the obtained results. Moreover, $\mathcal{P}$ (7) and
$\mathcal{Q}$ (7) turns into

$\mathcal{P}(z)=\frac{\frac{7}{2}\Omega_{m_{0}}(1+z)^{3}+4\Omega_{r_{0}}(1+z)^{4}+2\Omega_{\Lambda}}{(1+z)[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{\Lambda}]}$
$\mathcal{Q}(z)=\frac{3\Omega_{m_{0}}(1+z)+4\Omega_{r_{0}}(1+z)^{2}}{2[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{\Lambda}]}.$

Then, the GDE for null vector fields becomes

$\displaystyle\frac{d^{2}\eta}{dz^{2}}$
$\displaystyle+\frac{\frac{7}{2}\Omega_{m_{0}}(1+z)^{3}+4\Omega_{r_{0}}(1+z)^{4}+2\Omega_{\Lambda}}{(1+z)[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{\Lambda}]}\frac{d\eta}{dz}$
$\displaystyle+\frac{3\Omega_{m_{0}}(1+z)+4\Omega_{r_{0}}(1+z)^{2}}{2[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}+\Omega_{\Lambda}]}\eta=0\,.$

We set $\Omega_{\Lambda}=0$ and $\Omega_{m_{0}}+\Omega_{r_{0}}=1$ for the
original Mattig relation, so we have

$\displaystyle\frac{d^{2}\eta}{dz^{2}}$
$\displaystyle+\frac{\frac{7}{2}\Omega_{m_{0}}(1+z)^{3}+4\Omega_{r_{0}}(1+z)^{4}}{(1+z)[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}]}\frac{d\eta}{dz}$
$\displaystyle+\frac{3\Omega_{m_{0}}(1+z)+4\Omega_{r_{0}}(1+z)^{2}}{2[\Omega_{m_{0}}(1+z)^{3}+\Omega_{r_{0}}(1+z)^{4}]}\eta=0\,.$

This gives us a hint that the generalized Mattig relation in $f(Q)$-gravity
can be obtained from (33). In FLRW universe, the angular diametral distance
$D_{A}(z)$ is given by [15]

$D_{A}(z)=\sqrt{\left|\frac{dA(z)}{d\Omega}\right|}$

where $dA$ is the area of the object and $d\Omega$ is the solid angle. Thus,
from (33), the GDE in terms of the angular diametral distance is

$\frac{d^{2}D_{A}}{dz^{2}}+\mathcal{P}(H,\dot{H},z)\frac{dD_{A}}{dz}+\mathcal{Q}(H,\dot{H},z)D_{A}=0$
(36)

where $\mathcal{P}$ and $\mathcal{Q}$ is given in (7) and (7). This equation
satisfies the initial conditions (for $z\geq z_{0}$)

$\displaystyle D_{A}(z)|_{z=z_{0}}$ $\displaystyle=0$ (37)
$\displaystyle\frac{dD_{A}}{dz}(z)|_{z=z_{0}}$
$\displaystyle=\frac{H_{0}}{H(z_{0})(1+z_{0})}$ (38)

where $H(z_{0})$ is the modified Friedmann equation (31) at $z=z_{0}$.

##  8\. Dyer-Roeder like equation in $f(Q)$-gravity

Finally we describe a relation that is a tool to investigate cosmological
distances in inhomogeneous universes. The Dyer-Roeder equation is a
differential equation for the diametral angular distance as a function of the
redshift. The Dyer-Roeder equation in GR is given by [22]

$(1+z)^{2}\mathcal{F}(z)\frac{d^{2}D_{A}}{dz^{2}}+(1+z)\mathcal{G}(z)\frac{dD_{A}}{dz^{2}}+\mathcal{H}(z)D_{A}=0$

where

$\mathcal{F}(z)=H^{2}(z)\,,$
$\mathcal{G}(z)=(1+z)H(z)\frac{dH}{dz}+2H^{2}(z)\,,$

and

$\mathcal{H}(z)=\frac{3\tilde{\alpha}(z)}{2}\Omega_{m0}(1+z)^{3}\,,$

where $\tilde{\alpha}(z)$ is the smoothness parameter, which provides the
property of inhomogeneities in the energy density. The influence of the
smoothness parameter $\tilde{\alpha}$ in the behavior of $D_{A}(z)$ is
discussed in [15, 22]. Now, we express the Dyer-Roeder like equation in
$f(Q)$-gravity. First, notice that the terms involving the derivatives of
$D_{A}$ in (36) are from the transformation
$\frac{d}{d\nu}\rightarrow\frac{d}{dz}$, while the terms with $D_{A}$ are from
the Ricci focusing in (26). Then, define a mass-fraction $\tilde{\alpha}$ of
matter in the universe, and replacing the $\rho$ in the Ricci focusing with
$\tilde{\alpha}\rho$. Hence, from (33), and consider the case
$\Omega_{r_{0}}=0$, we obtain

$\displaystyle\frac{d^{2}D_{A}^{DR}}{dz^{2}}$
$\displaystyle+\frac{\frac{7}{2}\Omega_{m_{0}}(1+z)^{3}+2\Omega_{DE}-\frac{12\dot{H}f^{\prime\prime}}{f^{\prime}}[\Omega_{m_{0}}(1+z)^{3}+\Omega_{DE}]}{(1+z)[\Omega_{m_{0}}(1+z)^{3}+\Omega_{DE}]}\frac{dD_{A}^{DR}}{dz}$
$\displaystyle+\frac{3\tilde{\alpha}(z)\Omega_{m_{0}}(1+z)^{3}-\frac{24\dot{H}f^{\prime\prime}}{f^{\prime}}[\Omega_{m_{0}}(1+z)^{3}+\Omega_{DE}]}{2(1+z)^{2}[\Omega_{m_{0}}(1+z)^{3}+\Omega_{DE}]}dD_{A}^{DR}=0$
(39)

where $D_{A}^{DR}$ denote the Dyer-Roeder distance in $f(Q)$-gravity. This
equation also satisfies the conditions stated in (37) and (38). In the case of
$f(Q)=Q-2\Lambda$, (8) reduces to the standard form of GR.

##  9\. Conclusion

In the core of this paper lies the Ricci curvature tensor corresponding to the
Levi-Civita connection, expressed in terms of the tensor $Q_{\mu\nu\lambda}$
with the covariant form of the field equations of $f(Q)$-gravity theory. In
the FLRW universe, the GDE corresponding to these GR comparable terms of
$f(Q)$-gravity is acquired for the fundamental observer and the past-directed
null vector fields. The null vector field case provides an important results
which is the generalisation of the Mattig relation and the differential
equation for the diametral angular distance in $f(Q)$-gravity. As an
extension, the Dyer-Roeder equation was considered.

## References

  * [1] R. Akdrovandi and J. G. Pereira, Teleparallel Gravity: An Introduction. Springer-Verlag, Berlin, 2013.
  * [2] I. Ayuso, R. Lazkoz and V. Salzano, Observational constraints on cosmological solutions of f(Q) theories. Phys. Rev. D 103 (2021) 063505.
  * [3] P. Brax, What makes the Universe accelerate? A review on what dark energy could be and how to test it. Reports on Progress in Physics 81 (2018) 016902.
  * [4] D. L. Caceres, L. Castaneda and J. M. Tejeiro, Geodesic deviation equation in Bianchi cosmologies. J. Phys. Conf. Ser. 229 (2010) 012076.
  * [5] F. Darabi, M. Mousavi and K. Atazadeh, Geodesic deviation equation in f(T) gravity. Phys. Rev. D 91 (2015) 084023.
  * [6] G. F. R. Ellis and H. V. Elst, Deviation of geodesics in FLRW spacetime geometries. arXiv:gr-qc/9709060.
  * [7] A. De Felice and S. Tsujikawa, f(R) theories. Living Reviews Relativity 13 (2010) 3.
  * [8] A. Guarnizo, L. Castaneda, and J. M. Tejeiro, Erratum to: geodesic deviation equation in f(R) gravity. Gen. Rel. Grav. 47 (2015) 109.
  * [9] A. Guarnizo, L. Castaneda, and J. M. Tejeiro, Geodesic deviation equation in f(R) gravity. Gen. Rel. Grav. 43 (2011) 2713.
  * [10] K. Hayashi and T. Shirafuji, New General Relativity. Phys. Rev. D 19 (1979) 3524.
  * [11] L. Iorio and E. N. Saridakis, Einstein’s other gravity and the acceleration of the Universe. Mon. Not. R. Astron. Soc. 426 (2012) 1555.
  * [12] J. B. Jiménez, L. Heisenberg and T. S. Koivisto and S. Pekar, Cosmology in f(Q) geometry. Phys. Rev. D 101 (2020) 103507.
  * [13] J. M. Nester and H. J. Yo, Symmetric teleparallel general relativity. Chin. J. Phys. 37 (1999) 113.
  * [14] S. Nojiri, S. D. Odintsov and V. K. Oikonomou. Modified gravity theories on a nutshell: inflation, bounce and late-time evolution. Phys. Rept. 692 (2017) 1.
  * [15] P. Schneider, J. Ehlers and E. E. Falco, Gravitational lenses. Springer-Verlag, 1999.
  * [16] T. P. Sotiriou and V. Faraoni, f(R) theories of gravity. Rev. Mod. Physics. 82 (2010) 451.
  * [17] N. Tamanini and C. G. Boehmer, Good and bad tetrads in f(T) gravity. Phys. Rev. D 86 (2012) 044009.
  * [18] R. M. Wald, General Relativity. University of Chicago Press, Chicago, 1984.
  * [19] C. M. Will, The confrontation between General Relativity and experiment. Living Reviews Relativity 17 (2014) 4.
  * [20] Y. Xu, G. Li, T. Harko and S. D. Liang, f(Q,T) gravity. Eur. Phys. J. C 70 (2019) 708.
  * [21] D. Zhao, Covariant formulation of f(Q) theory. arXiv:2104.02483[gr-qc]
  * [22] T. Okamura and T. Futamase, Distance-Redshift Relation in a Realistic Inhomogeneous Universe. Prog. Theor. Phys. 122 (2019) 511.

##  10\. Appendix

From (2), we can get all the non-metricity tensors as follow

$\displaystyle Q_{\lambda\mu\nu}$ $\displaystyle=\nabla_{\lambda}g_{\mu\nu}$
$\displaystyle Q^{\lambda}{}_{\mu\nu}$
$\displaystyle=g^{\lambda\rho}Q_{\rho\mu\nu}=g^{\lambda\rho}\nabla_{\rho}g_{\mu\nu}=\nabla^{\lambda}g_{\mu\nu}$
$\displaystyle Q_{\lambda}{}^{\mu}{}_{\nu}$
$\displaystyle=g^{\mu\rho}Q_{\lambda\rho\nu}=g^{\mu\rho}\nabla_{\lambda}g_{\rho\nu}=-g_{\rho\nu}\nabla_{\lambda}g^{\mu\rho}$
$\displaystyle Q_{\lambda\mu}{}^{\nu}$
$\displaystyle=g^{\nu\rho}Q_{\lambda\mu\rho}=g^{\nu\rho}\nabla_{\lambda}g_{\mu\rho}=-g_{\mu\rho}\nabla_{\lambda}g^{\nu\rho}$
$\displaystyle Q^{\lambda\mu}{}_{\nu}$
$\displaystyle=g^{\lambda\rho}g^{\mu\gamma}\nabla_{\rho}g_{\gamma\nu}=g^{\mu\gamma}\nabla^{\lambda}g_{\gamma\nu}=-g_{\gamma\nu}\nabla^{\lambda}g^{\mu\gamma}$
$\displaystyle Q^{\lambda}{}_{\mu}{}^{\nu}$
$\displaystyle=g^{\lambda\rho}g^{\nu\gamma}\nabla_{\rho}g_{\mu\gamma}=g^{\nu\gamma}\nabla^{\lambda}g_{\mu\gamma}=-g_{\mu\gamma}\nabla^{\lambda}g^{\nu\gamma}$
$\displaystyle Q_{\lambda}{}^{\mu\nu}$
$\displaystyle=g^{\mu\rho}g^{\nu\gamma}\nabla_{\lambda}g_{\rho\gamma}=-g^{\mu\rho}g_{\rho\gamma}\nabla_{\lambda}g^{\nu\gamma}=-\nabla_{\lambda}g^{\mu\nu}$
$\displaystyle Q^{\lambda\mu\nu}$
$\displaystyle=-\nabla^{\lambda}g^{\mu\nu}\,.$

Recall in (7), we have

$Q=-\frac{1}{4}Q_{\lambda\mu\nu}Q^{\lambda\mu\nu}+\frac{1}{2}Q_{\lambda\mu\nu}Q^{\mu\lambda\nu}+\frac{1}{4}Q_{\lambda}Q^{\lambda}-\frac{1}{2}Q_{\lambda}\tilde{Q}^{\lambda}\,.$

By using the above results and the FLRW metric in (17), we obtain

$\displaystyle Q_{\lambda\mu\nu}Q^{\lambda\mu\nu}$
$\displaystyle=-\nabla_{\lambda}g_{\mu\nu}\nabla^{\lambda}g^{\mu\nu}=-12H^{2}$
$\displaystyle Q_{\lambda\mu\nu}Q^{\mu\lambda\nu}$
$\displaystyle=-\nabla_{\lambda}g_{\mu\nu}\nabla^{\mu}g^{\lambda\nu}=0$
$\displaystyle Q_{\lambda}Q^{\lambda}$
$\displaystyle=(g_{\mu\rho}\nabla_{\lambda}g^{\mu\rho})(g_{\nu\gamma}\nabla^{\lambda}g^{\nu\gamma})=-36H^{2}$
$\displaystyle Q_{\lambda}\tilde{Q}^{\lambda}$
$\displaystyle=(g_{\mu\rho}\nabla_{\lambda}g^{\mu\rho})(\nabla_{\nu}g^{\lambda\nu})=0\,.$

Hence, we have

$Q=-\frac{1}{4}(-12H^{2})+\frac{1}{4}(-36H^{2})=-6H^{2}\,.$

Next, we try to simplify the terms

$\frac{1}{2f^{\prime}}\bigg{[}(\delta^{\alpha}_{\gamma}\mathcal{D}_{\delta\beta}-\delta^{\alpha}_{\delta}\mathcal{D}_{\gamma\beta}+g_{\beta\delta}\mathcal{D}^{\alpha}_{\gamma}-g_{\beta\gamma}\mathcal{D}^{\alpha}_{\delta})f^{\prime}V^{\beta}V^{\delta}\bigg{]}\eta^{\gamma}$
$\frac{4}{3}\epsilon\eta^{\alpha}f^{\prime\prime}P^{\lambda\nu}{}_{\nu}\nabla_{\lambda}Q$

as stated in (4). From the definition of the operator $\mathcal{D}_{\mu\nu}$
in (14), we have

$\displaystyle\delta^{\alpha}_{\gamma}V^{\beta}V^{\delta}\eta^{\gamma}\mathcal{D}_{\delta\beta}f^{\prime}$
$\displaystyle=-2\eta^{\alpha}V^{\delta}V^{\beta}P^{\lambda}{}_{\delta\beta}\nabla_{\lambda}Qf^{\prime\prime}$
$\displaystyle-\delta^{\alpha}_{\delta}V^{\beta}V^{\delta}\eta^{\gamma}\mathcal{D}_{\gamma\beta}f^{\prime}$
$\displaystyle=2V^{\alpha}\eta^{\gamma}V^{\beta}P^{\lambda}{}_{\gamma\beta}\nabla_{\lambda}Qf^{\prime\prime}$
$\displaystyle
g_{\beta\delta}V^{\beta}V^{\delta}\eta^{\gamma}\mathcal{D}^{\alpha}_{\gamma}f^{\prime}$
$\displaystyle=-2\epsilon\eta^{\gamma}P^{\lambda\alpha}{}_{\gamma}\nabla_{\lambda}Qf^{\prime\prime}$
$\displaystyle-
g_{\beta\gamma}V^{\beta}V^{\delta}\eta^{\gamma}\mathcal{D}^{\alpha}_{\delta}f^{\prime}$
$\displaystyle=0\,.$

Since $Q$ is only time-dependent, so the summation of the index $\lambda$
reduces to only the 0 component. Then, we verify the terms as follow

$\displaystyle P^{0}{}_{\mu\nu}$
$\displaystyle=\frac{1}{4}\left(-Q^{0}{}_{\mu\nu}+Q_{\mu}{}^{0}{}_{\nu}+Q_{\nu}{}^{0}{}_{\mu}+Q^{0}g_{\mu\nu}-\tilde{Q}^{0}g_{\mu\nu}-\frac{1}{2}\delta^{0}_{\mu}Q_{\nu}-\frac{1}{2}\delta^{0}_{\nu}Q_{\mu}\right)$
$\displaystyle=\frac{1}{4}\left(\nabla_{0}g_{\mu\nu}-6Hg_{\mu\nu}+\frac{1}{2}\delta^{0}_{\mu}g_{\alpha\beta}\nabla_{\nu}g^{\alpha\beta}+\frac{1}{2}\delta^{0}_{\nu}g_{\alpha\beta}\nabla_{\mu}g^{\alpha\beta}\right)$
$\displaystyle P^{0\mu}{}_{\nu}$
$\displaystyle=\frac{1}{4}\left(-Q^{0\mu}{}_{\nu}+Q^{\mu
0}{}_{\nu}+Q_{\nu}{}^{0\mu}+Q^{0}\delta^{\mu}_{\nu}-\tilde{Q}^{0}\delta^{\mu}_{\nu}-\frac{1}{2}g^{0\mu}Q_{\nu}-\frac{1}{2}\delta^{0}_{\nu}Q^{\mu}\right)$
$\displaystyle=\frac{1}{4}\left(-g_{\rho\nu}\nabla_{0}g^{\mu\rho}-6H\delta^{\mu}_{\nu}+\frac{1}{2}g^{0\mu}g_{\alpha\beta}\nabla_{0}g^{\alpha\beta}+\frac{1}{2}\delta^{0}_{\nu}g_{\rho\nu}g^{\alpha\mu}\nabla_{\alpha}g^{\rho\nu}\right)$
$\displaystyle P^{0\nu}{}_{\nu}$
$\displaystyle=\frac{1}{4}\left(-g_{\rho\nu}\nabla_{0}g^{\rho\nu}-6H\delta^{\nu}_{\nu}+\frac{1}{2}g^{0\nu}g_{\alpha\beta}\nabla_{0}g^{\alpha\beta}+\frac{1}{2}g^{0\nu}g_{\rho\nu}\nabla_{\nu}g^{\rho\nu}\right)$
$\displaystyle=\frac{1}{4}(2g_{\alpha\beta}\nabla_{0}g^{\alpha\beta})$
$\displaystyle=-3H\,.$

Notice that if $\mu\neq\nu$, then $P^{0}{}_{\mu\nu}=0$. This implies that

$\displaystyle V^{\mu}V^{\nu}P^{0}{}_{\mu\nu}$
$\displaystyle=V^{0}V^{0}P^{0}{}_{00}+\sum_{k}V^{k}V^{k}P^{0}{}_{kk}$
$\displaystyle=\frac{1}{4}\left[V^{0}V^{0}(\nabla_{0}g_{00}-6Hg_{00}-6H)+\sum_{k}V^{k}V^{k}(\nabla_{0}g_{kk}-6Hg_{kk})\right]$
$\displaystyle=\frac{1}{4}\left[V^{0}V^{0}(0)+\sum_{k}V^{k}V^{k}\nabla_{0}g_{kk}-(6H)\sum_{k}V^{k}V^{k}g_{kk}\right]$
$\displaystyle=\frac{1}{4}\left[\sum_{k}V^{k}V_{k}g^{kk}\nabla_{0}g_{kk}-(6H)\sum_{k}V^{k}V_{k}\right]$
$\displaystyle=\frac{1}{4}\left[\sum_{k}V^{k}V_{k}(2H-6H)\right]$
$\displaystyle=-\frac{1}{4}(\epsilon-V_{0}V^{0})4H$
$\displaystyle=-H(\epsilon+E^{2})$
$\displaystyle\eta^{\mu}V^{\nu}P^{0}{}_{\mu\nu}$
$\displaystyle=\eta^{0}V^{0}P^{0}{}_{00}+\sum_{k}\eta^{k}V^{k}P^{0}{}_{kk}$
$\displaystyle=\frac{1}{4}\sum_{k}\eta^{k}V^{k}(\nabla_{0}g_{kk}-6Hg_{kk})$
$\displaystyle=\frac{1}{4}\sum_{k}\eta^{k}V_{k}(2H-6H)$ $\displaystyle=0$
$\displaystyle\eta^{\nu}P^{0\mu}{}_{\nu}$
$\displaystyle=\eta^{0}P^{0\mu}{}_{0}+\eta^{i}P^{0\mu}{}_{i}$
$\displaystyle=\frac{1}{4}\eta^{i}(-g_{\rho
i}\nabla_{0}g^{\rho\mu}-\delta^{\mu}_{i}6H)$
$\displaystyle=\frac{1}{4}(-\eta_{\rho}\nabla_{0}g^{\rho\mu}-\eta^{\mu}6H)$
$\displaystyle=\frac{1}{4}(2\eta^{\mu}H-\eta^{\mu}H)$
$\displaystyle=-H\eta^{\mu}\,.$

Thus, we have

$\displaystyle-2\eta^{\alpha}V^{\delta}V^{\beta}P^{\lambda}{}_{\delta\beta}\nabla_{\lambda}Qf^{\prime\prime}$
$\displaystyle=-2\eta^{\alpha}(-H)(\epsilon+E^{2})(-12H\dot{H})f^{\prime\prime}$
$\displaystyle=-24H^{2}\dot{H}f^{\prime\prime}(\epsilon+E^{2})\eta^{\alpha}$
$\displaystyle
2V^{\alpha}\eta^{\gamma}V^{\beta}P^{\lambda}{}_{\gamma\beta}\nabla_{\lambda}Qf^{\prime\prime}$
$\displaystyle=0$
$\displaystyle-2\epsilon\eta^{\gamma}P^{\lambda\alpha}{}_{\gamma}\nabla_{\lambda}Qf^{\prime\prime}$
$\displaystyle=-2\epsilon(-H\eta^{\alpha})(-12H\dot{H})f^{\prime\prime}$
$\displaystyle=-24H^{2}\dot{H}f^{\prime\prime}\epsilon\eta^{\alpha}\,.$

Therefore,

$\frac{1}{2f^{\prime}}\bigg{[}(\delta^{\alpha}_{\gamma}\mathcal{D}_{\delta\beta}-\delta^{\alpha}_{\delta}\mathcal{D}_{\gamma\beta}+g_{\beta\delta}\mathcal{D}^{\alpha}_{\gamma}-g_{\beta\gamma}\mathcal{D}^{\alpha}_{\delta})f^{\prime}V^{\beta}V^{\delta}\bigg{]}\eta^{\gamma}=\frac{1}{2f^{\prime}}[-24H^{2}\dot{H}f^{\prime\prime}(2\epsilon+E^{2})]\eta^{\alpha}$

and

$\displaystyle\frac{4}{3}\eta^{\alpha}f^{\prime\prime}P^{\lambda\nu}{}_{\nu}\nabla_{\lambda}Q$
$\displaystyle=\frac{4}{3}\epsilon\eta^{\alpha}f^{\prime\prime}(-3H)(-12H\dot{H})$
$\displaystyle=48H^{2}\dot{H}f^{\prime\prime}\epsilon\eta^{\alpha}\,.$


|  | HSF-DOC-2022-01
---|---|---
|  | May 17, 2022
|  | Copyright (C) 2022 CERN, Princeton and Fermilab, licence CC-BY-4.0

# The HEP Software Foundation Community

The HEP Software Foundation
Contact editors:
Graeme A Stewart, CERN<EMAIL_ADDRESS>
Peter Elmer, Princeton University<EMAIL_ADDRESS>
Elizabeth Sexton-Kennedy, Fermilab<EMAIL_ADDRESS>
(April 2022)

###### Abstract

The HEP Software Foundation was founded in 2014 to tackle common problems of
software development and sustainability for high-energy physics. In this paper
we outline the motivation for the founding of the organisation and give a
brief history of its development. We describe how the organisation functions
today and what challenges remain to be faced in the future.

##  1 History

Over the past 50 years, the experimental particle, nuclear and astroparticle
physics communities have iteratively evolved a significant amount of community
structure. This was a natural result of the growing size, scale and time
duration of the experiments and the centralisation of facilities at large
laboratories. National, and now international, collaborations are typically
required to build, operate and maintain the large detectors used in these
experiments. No single university or laboratory can provide all of the
necessary expertise and required effort. The largest collaborations have grown
from 100s of collaborators in the 1990s to 1000s at (for example) the Large
Hadron Collider (LHC) at CERN. This community has also developed a broad
ecosystem of methodologies and technologies that are used to build successive
experiments and upgrades to existing experiments. While a specific instrument
can necessarily only be used for a single experiment at any given time, the
large commonalities in methodology should permit the development of research
software which can be used widely in the community for different projects.
Despite this, much of the software development remained somewhat siloed within
individual experiments or, at most, one or another host laboratory, with only
a few exceptions.

This is not to say that in the history of HEP there was no common software.
CERNLIB CERNLIB (1) was a foundation library written in the Fortran era and
used by many experiments. Elements of it have been rewritten in C++ and
constitute some of the most widely used software packages in the field.
Projects such as ROOT, Geant4, and various generators have effectively acted
as common glue for many experiments. However in the software layers above
these foundation and toolkit libraries, redundant solutions that are difficult
to evolve and sustain over time (years or decades for large experiments!) are
common. To a large extent, software speed, performance and efficiency had been
ignored previously, because the costs due to inefficient software could be
ignored in the past. Software is as much an intellectual product as well as a
tool, thus a new approach was needed.

First steps in the direction of collaborating on modernising, in 2011-2012,
led to the formation of a cross-experiment “Concurrency Forum” CERN-RD-
MULTICORE (2) to discuss the specific software challenges brought by changes
in microprocessor technology (multi-core, wide vectors, GPUs). Driven
initially by CERN and Fermilab, the forum demonstrated community interest in
wider software collaborations. By 2014-2015, a number of colleagues involved
in HEP software for many years were discussing a more ambitious and broader
scope for research software collaborations in HEP. This eventually led to the
formation of the High-Energy Physics (HEP) Software Foundation (HSF). The
driving motivations for this initiative were that the physics upgrades
anticipated in the coming decades, particularly the High-Luminosity LHC HL-LHC
(3), would put enormous pressure on the software used in HEP; that much of our
software was already decades old; the changes in microprocessor technology
brought new challenges to the table; and that there was an urgent need to
train new talent and to attract investment to the field, which could be better
supported when common, multi-experiment, efforts were promoted. More
generally, additional community structure which promotes research software
collaborations, not tied to single experiments or laboratories, has greater
potential to enhance the longer term sustainability of the software.

The very first workshop CERN-WS (4) attempted to build on the community
experience within the large experiments, however there was too much discussion
of “governance” questions. Individual experiments need to operate a large
well-integrated detector, manage pooled resources (such as computing and
storage) and at the end of the day produce scientific publications signed by
the entire collaboration. Thus governance questions within experiments are
critical. This top-down approach was less obvious for the envisioned research
software collaborations, which can be more “ecosystem-like”. It also made
engaging experiments of very different sizes more challenging. By the end of
the workshop most participants had concluded that a different structure was
needed. Subsequent workshops, one in North America and one in Europe to aid
inclusivity SLAC-WS (5, 6), brought together many HEP experiments, HEP
specific software projects and other non-HEP organisations, such as the Apache
Software Foundation and the Software Sustainability Institute. Here, the
principle of a “do-ocracy” was enshrined, to encourage activity from the
bottom-up, with developers leading, which was a far more productive approach
to community building.

From these early workshops the idea was born of preparing a Community White
Paper (CWP), laying out the roadmap for software and computing in HEP in the
2020s. As a way to fortify these community-led efforts, the Worldwide LHC
Computing Grid (WLCG) gave a formal charge to the HSF to prepare this paper
CWP-Charge (7). Many individuals volunteered and the US National Science
Foundation was an early investor with dedicated funding to help organise and
carry out CWP workshops SDSC-CWP (8, 9). These workshops were focal points of
an intense year of identifying key domain areas and potential solutions for
the field and setting up working groups who would prepare topic-specific white
papers. This approach was able to greatly broaden the involvement of key
individuals and solidified the active communities around the HSF. Once the
working groups had produced their domain-specific white papers, an editorial
team took charge of synthesising a unified single version of the white paper,
which was published with extremely wide community support, 310 signing authors
from 124 institutes Albrecht2019 (10).

The CWP was not the only activity happening in the HSF at this time.
Particularly where there was an identified gap between different communities,
be these inter-experiment, between projects, or even between theory and
experiment, the HSF was ideally placed to bridge gaps and bring different
people together. Workshops on analysis ecosystems ANALYSIS-ECO (11) and on
computational aspects of event generators COMP-GEN (12) were notable examples.
In addition the HSF was a natural forum for considering more radical software
developments and their potential, gathering experiment and expert feedback on
progress and possibilities GEANTV-RD (13).

##  2 HSF Activities

While much of the early HSF activity was focused on software community
building through the CWP, working groups were started where common topics were
obviously identified within the community, e.g., in the domain of packaging
and distributing HEP software stacks. These groups brought together experts
and interested parties and focused around technical discussions.

In the wake of the CWP it was clear that this model would work very well for
the areas of most concern for the future, so the model was broadened and, over
a few years, eight working groups were established with about 3 conveners in
each case, appointed annually and with a nomination system that allows both
stakeholder (experiment, institution) input and bottom-up volunteers for
running these groups. In order to marshall these activities it was necessary
to have some critical amount of binding effort, so the decision of CERN
management to allow significant (0.5 FTE) time from one person to the HSF was
crucial and has had a key multiplication effect.

Where these groups are involved in areas that are ‘traditional’, e.g.,
detector simulation or reconstruction, there is a strong involvement with
ongoing developments in the experiments and in well established software
projects. The focus is the exchange of ideas and discussions of common
problems. In a number of cases the HSF had identified topics of interest that
were simply not covered elsewhere in the field and then the HSF working group
has had a further leading and catalysing effect. This is particularly the case
for the use of Python in HEP, led by the PyHEP working group; and for
computational aspects of physics event generators, led by the Generators
working group Valassi2021 (14).

In addition to being a focus for the exchange of ideas and techniques, when
WGs identify a concrete topic where a paper can usefully be prepared, the HSF
is a natural place for organising pan-experiment input and encouraging some
standardisation (e.g., in analysis level metadata or in detector conditions
data access). This has been recognised by more formal bodies outside the HSF,
such as WLCG and the LHCC, who often ask the HSF to marshall community inputs
for updates on development activities and reviews of development plans. This
is an important point in terms of recognising the contribution of the
organisation. In turn, this fosters recognition for the contribution of our
individual members and helps their careers to advance.

This engagement with strategic bodies in the field, where the HSF advocates
for software investment stewart_graeme_a_2018_2413005 (15), leads the HSF to
work in tandem with other funded software R&D projects in HEP, where the HSF
will support funding applications and then work with these projects, enhancing
their connection to the community and the impact of their work. Practically
projects can contribute through the working groups closest to their R&D areas.

The HSF also engages with other like minded bodies and collaborates on regular
series of meetings regarding accelerator programming or software and computing
for nuclear physics and also organises itself as an umbrella organisation,
with CERN, to run Google Summer of Code for HEP software projects and
experiments.

In the past the HSF organised face-to-face workshops and the intention is to
restart such activities once the pandemic passes, but in the meantime virtual
workshops have played a role. In some circumstances these can even have a
greater impact, with the PyHEP workshops in 2020 and 2021 registering more
than 1000 participants PYHEP20 (16, 17). In large part this reflects a strong
didactic element in the Python area. This thread is reflected also in that the
HSF, together with IRIS-HEP, SIDIS, The Carpentries and the ROOT project IRIS-
HEP (18, 19, 20, 21, 22), has put a strong emphasis on training activities
Malik2021 (23) and now runs regular training events in fundamental software
skills and in C++ programming. This is seen as a critical activity in the
field and attempts are now also being made to have these activities as feeders
to encourage trainees to be involved in software development and training.

##  3 Outcomes and Conclusions

Contrary to initial ideas, the HSF has not, by and large, run software
projects themselves - without actually having resources to disburse it is
better to allow such projects to be independent and work with the HSF as is
useful. That said, HSF events have proved to be fertile ground for people from
different backgrounds to meet (e.g., nuclear and astroparticle physics) and
even to start common software projects that then take on a life of their own.
Fostering funded projects has led to new investment in software R&D in HEP and
a higher recognition of the importance of promoting excellent software
developers in their careers; this is a major success that was a direct outcome
of the CWP process.

The HSF has now established itself as a recognised part of the HEP software
landscape where it links strategic bodies to the community of software
developers. It remains a challenge to continue to build the next generation of
HEP software developers and make them feel involved and part of an
organisation like the HSF, but work to improve training and engage younger
colleagues through the working group process is hoped to improve this. Looking
forward to post-pandemic activities, where face-to-face interactions can
happen again, will also help to continue to build HEP software communities.

## References

  * (1) “CERN Program Library”  URL: https://en.wikipedia.org/wiki/CERN_Program_Library
  * (2) “Forum on Concurrent Programming Models and Frameworks”  URL: https://concurrency.web.cern.ch/concurrency/index.html
  * (3) “The High-Luminosity LHC project”  URL: https://home.cern/science/accelerators/high-luminosity-lhc
  * (4) “HEP Software Collaboration meeting”, 2014  URL: https://indico.cern.ch/event/297652/
  * (5) “HEP Software Foundation Workshop (SLAC)”, 2015  URL: https://indico.cern.ch/event/357737/
  * (6) “HEP Software Foundation Workshop (LAL)”, 2016  URL: https://indico.cern.ch/event/496146/
  * (7) “Charge for Producing a HSF Community White Paper”, 2016  URL: https://hepsoftwarefoundation.org/assets/CWP-Charge-HSF.pdf
  * (8) “HEP Software Foundation Workshop (SDSC)”, 2017  URL: https://indico.cern.ch/event/570249/
  * (9) “HEP Software Foundation Workshop (LAPP)”, 2017  URL: https://indico.cern.ch/event/613093/
  * (10) Johannes Albrecht et al.  “A Roadmap for HEP Software and Computing R&D for the 2020s”  In _Computing and Software for Big Science_ 3.1, 2019, pp. 7  DOI: 10.1007/s41781-018-0018-8
  * (11) “HEP Analysis Ecosystem Workshop”  URL: https://indico.cern.ch/event/570249/
  * (12) “Physics Event Generator Computing Workshop”  URL: https://indico.cern.ch/event/751693/
  * (13) “HEP Software Community Meeting on GeantV R&D”  URL: https://indico.cern.ch/event/570876/
  * (14) Andrea Valassi et al.  “Challenges in Monte Carlo Event Generator Software for High-Luminosity LHC”  In _Computing and Software for Big Science_ 5.1, 2021, pp. 12  DOI: 10.1007/s41781-021-00055-1
  * (15) Graeme A Stewart  “The Importance of Software and Computing to Particle Physics”  Zenodo, 2018  DOI: 10.5281/zenodo.2413005
  * (16) “PyHEP 2020 (virtual) Workshop”, 2020  URL: https://indico.cern.ch/e/pyhep2020
  * (17) “PyHEP 2021 (virtual) Workshop”, 2021  URL: https://indico.cern.ch/e/pyhep2021
  * (18) “Institute for Research and Innovation in Software for High Energy Physics (IRIS-HEP)”  URL: https://iris-hep.org/
  * (19) “Software Institute for Data-Intensive Sciences”  URL: https://sidis.web.cern.ch/
  * (20) “The Carpentries”  URL: https://carpentries.org/
  * (21) R. Brun and F. Rademakers  “ROOT: An object oriented data analysis framework”  In _New computing techniques in physics research V. Proceedings, 5th International Workshop, AIHENP ’96, Lausanne, Switzerland, September 2-6, 1996_ A389, 1997, pp. 81–86  DOI: 10.1016/S0168-9002(97)00048-X
  * (22) “ROOT Data Analysis Framework”  URL: https://root.cern/
  * (23) Sudhir Malik et al.  “Software Training in HEP”  In _Computing and Software for Big Science_ 5.1, 2021, pp. 22  DOI: 10.1007/s41781-021-00069-9


$\displaystyle h(x)=100(x_{1}^{2}-x_{2})^{2}+(x_{3}-1)^{2}+(x_{1}-1)^{2}.$
$\displaystyle-10\leq x_{i}\leq 10,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \forall i=1,\ldots,4.$

The global maximum is $0$ with $x^{\star}=[1,1,1,1]^{\top}$.

###  C.5 Friedman

The Friedman problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-\left(10y_{1}+20(x_{3}-0.5)^{2}+10x_{4}+5x_{5}\right).$
$\displaystyle h(x)=\sin(\pi x_{1}x_{2}).$ $\displaystyle 0\leq x_{i}\leq
1,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall
i=1,\ldots,5.$

The global maximum is $27.5$ with
$x^{\star}=[x_{1}^{\star},x_{2}^{\star},0.5,-1.5,-1.5]^{\top}$ for any
$x_{1}^{\star}$ and $x_{2}^{\star}$ satisfying $\sin(\pi
x_{1}^{\star}x_{2}^{\star})=-1$.

###  C.6 Dolan

The Dolan problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-\left(y_{1}-y_{2}+0.2x_{5}^{2}-x_{2}-1\right).$ $\displaystyle
h(x)=\begin{bmatrix}(x_{1}+1.7x_{2})\sin(x_{1})\\\
1.5x_{3}-0.1x_{4}\cos(x_{5}+x_{4}-x_{1})\end{bmatrix}.$ $\displaystyle-100\leq
x_{i}\leq 100,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\forall i=1,\ldots,5.$

The global maximum is $529.87$ with
$x^{\star}=[98.964,100,100,99.224,-0.25]^{\top}$.

###  C.7 Rosenbrock

The Rosenbrock problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-\bigg{(}\textstyle\sum_{i=1}^{3}(100y_{i}^{2}+(1-x_{i})^{2})+100(x_{5}-x_{4}^{2})+y_{4}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+100(x_{6}-x_{5}^{2})+(1-x_{5})^{2}\bigg{)}.$ $\displaystyle
h(x)=\begin{bmatrix}x_{2}^{2}-x_{1}^{2}\\\ x_{3}^{2}-x_{2}^{2}\\\
x_{4}^{2}-x_{3}^{2}\\\ (1-x_{4})^{2}\end{bmatrix}.$ $\displaystyle-2\leq
x_{i}\leq 2,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\forall i=1,\ldots,6.$

The global maximum is 0 with $x^{\star}=[0,0,0,0,0,0]^{\top}$.

###  C.8 Zakharov

The Zakharov problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-\left(\textstyle\sum_{i=}^{7}x_{i}^{2}+\textstyle\sum_{i=}^{7}(0.5ix_{i})^{2}+y_{1}\textstyle\sum_{i=}^{7}(0.5ix_{i})^{2}\right).$
$\displaystyle h(x)=\textstyle\sum_{i=}^{7}(0.5ix_{i})^{2}.$
$\displaystyle-5\leq x_{i}\leq 10,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \forall i=1,\ldots,7.$

The global maximum is $0$ with $x^{\star}=[0,0,0,0,0,0,0]^{\top}$.

###  C.9 Powell

The Powell problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-\bigg{(}y_{1}+(x_{5}+10x_{6})^{2}+y_{2}+5(x_{7}-x_{8})^{2}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+(x_{2}-2x_{3})^{4}+y_{3}+10(x_{1}-x_{4})^{4}+y_{4}\bigg{)}.$ $\displaystyle
h(x)=\begin{bmatrix}(x_{1}+10x_{2})^{2}\\\ 5(x_{3}-x_{4})^{2}\\\
(x_{6}-2x_{7})^{4}\\\ 10(x_{5}-x_{8})^{4}\end{bmatrix}.$ $\displaystyle-4\leq
x_{i}\leq 5,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\forall i=1,\ldots,8.$

The global maximum is $0$ with $x^{\star}=[0,0,0,0,0,0,0,0]^{\top}$.

###  C.10 Styblinski-Tang

The Styblinski-Tang problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-\left(\textstyle\sum_{i=1}^{4}y_{i}+\textstyle\sum_{i=5}^{9}(0.5x_{i}^{4}-16x_{i}^{2}+5x_{i})\right).$
$\displaystyle h(x)=\begin{bmatrix}0.5(x_{1}^{4}-16x_{1}^{2}+5x_{1})\\\
0.5(x_{2}^{4}-16x_{2}^{2}+5x_{2})\\\ 0.5(x_{3}^{4}-16x_{3}^{2}+5x_{3})\\\
0.5(x_{4}^{4}-16x_{4}^{2}+5x_{4})\end{bmatrix}.$ $\displaystyle-5\leq
x_{i}\leq 5,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\forall i=1,\ldots,9.$

The global maximum is $352.49$ with
$x^{\star}=-2.904[1,1,1,1,1,1,1,1,1]^{\top}$.

##  Appendix D Appendix: Constrained Synthetic Test Problems

###  D.1 Bazaraa

The Bazaraa problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-\left(2x_{1}^{2}+2x_{2}^{2}-y_{2}\right),$ $\displaystyle
g_{1}(x,y)$ $\displaystyle=-\left(5x_{1}+x_{2}-5\right),$ $\displaystyle
g_{2}(x,y)$ $\displaystyle=-\left(y_{1}-x_{1}\right).$ $\displaystyle
h(x)=\begin{bmatrix}2x_{2}^{2}\\\ 2x_{1}x_{2}+6x_{1}+4x_{2}\end{bmatrix}.$
$\displaystyle 0.01\leq x_{i}\leq 1,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \forall i=1,2.$

The global maximum is $6.613$ with $x^{\star}=[0.868,0.659]^{\top}$.

###  D.2 Spring

The Spring problem is defined by the following functions

$\displaystyle g_{0}(x,y)$ $\displaystyle=-\left(2y_{1}+2x_{3}\right),$
$\displaystyle g_{1}(x,y)$ $\displaystyle=-\left(2x_{2}^{2}-x_{1}\right),$
$\displaystyle g_{2}(x,y)$
$\displaystyle=-\left(\frac{4x_{2}^{2}-x_{1}x_{2}}{12566(x_{1}^{3}x_{2}-x_{1}^{4})}+\frac{1}{5108x_{1}x_{2}-1}\right),$
$\displaystyle g_{3}(x,y)$
$\displaystyle=-\left(1-140.45\frac{x_{1}}{y_{2}}\right),$ $\displaystyle
g_{4}(x,y)$ $\displaystyle=-\left(\frac{2}{3}(x_{1}+x_{2})-1\right).$
$\displaystyle h(x)=\begin{bmatrix}x_{1}^{2}x_{2}\\\
x_{2}^{3}x_{3}\end{bmatrix}.$ $\displaystyle 0.05$ $\displaystyle\leq
x_{1}\leq 2,$ $\displaystyle 0.25$ $\displaystyle\leq x_{2}\leq 1.3$
$\displaystyle 2$ $\displaystyle\leq x_{3}\leq 15.$

The global maximum is $-0.0127$ with $x^{\star}=[0.052,0.357,11.289]^{\top}$.

###  D.3 Ex314

The Ex314 problem is defined by the following functions

$\displaystyle g_{0}(x,y)$ $\displaystyle=y_{2},$ $\displaystyle g_{1}(x,y)$
$\displaystyle=-x_{1}y_{1}-2x_{2}^{2}+2x_{1}x_{2}+2x_{2}x_{3}-2x_{1}x_{3}-2x_{3}^{2}+20x_{1}-9x_{2}+13x_{3}-24,$
$\displaystyle g_{2}(x,y)$ $\displaystyle=x_{1}+x_{3}+x_{3}-4,$ $\displaystyle
g_{3}(x,y)$ $\displaystyle=3x_{2}+x_{6}-6.$ $\displaystyle
h(x)=\begin{bmatrix}4x_{1}-2x_{2}+2x_{3}\\\ x_{2}-x_{3}-2x_{1}\end{bmatrix}.$
$\displaystyle-2$ $\displaystyle\leq x_{1}\leq 2,$ $\displaystyle 0$
$\displaystyle\leq x_{2}\leq 6$ $\displaystyle-3$ $\displaystyle\leq x_{3}\leq
3.$

The global maximum is $4$ with $x^{\star}=[0.5,0.0,3.0]^{\top}$.

###  D.4 Rosen-Suzuki

The Rosen-Suzuki problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-(x_{1}^{2}+x_{2}^{2}+x_{4}^{2}-5x_{1}-5x_{2}+y_{1}),$
$\displaystyle g_{1}(x,y)$
$\displaystyle=8-x_{1}^{2}-x_{2}^{2}-x_{3}^{2}-x_{4}^{2}-x_{1}+x_{2}-x_{3}+x_{4},$
$\displaystyle g_{2}(x,y)$
$\displaystyle=10-x_{1}^{2}-2x_{2}^{2}-y_{2}+x_{1}+x_{4},$ $\displaystyle
g_{3}(x,y)$
$\displaystyle=5-2x_{1}^{2}-x_{2}^{2}-x_{3}^{2}-2x_{1}+x_{2}+x_{4}.$
$\displaystyle h(x)=\begin{bmatrix}2x_{3}^{2}-21x_{3}+7x_{4}\\\
x_{3}^{2}+2x_{4}^{2}\end{bmatrix}.$ $\displaystyle-2\leq x_{i}\leq
2,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall
i=1,\ldots,4.$

The global maximum is $44$ with $x^{\star}=[0,1,2,-1]^{\top}$.

###  D.5 st_bpv1

The st_bpv1 problem is defined by the following functions

$\displaystyle g_{0}(x,y)$ $\displaystyle=-(y_{1}+x_{2}x_{4}),$ $\displaystyle
g_{1}(x,y)$ $\displaystyle=-(30-y_{2}),$ $\displaystyle g_{2}(x,y)$
$\displaystyle=-(20-y_{3}),$ $\displaystyle g_{3}(x,y)$
$\displaystyle=-(x_{3}+x_{4}-15).$ $\displaystyle
h(x)=\begin{bmatrix}x_{1}x_{3}\\\ x_{1}+3x_{2}\\\ 2x_{1}+x_{2}\end{bmatrix}.$
$\displaystyle 0\leq x_{1}\leq 27,$ $\displaystyle 0\leq x_{2}\leq 16,$
$\displaystyle 0\leq x_{3}\leq 10,$ $\displaystyle 0\leq x_{4}\leq 10.$

The global maximum is $-10$ with $x^{\star}=[27,1,0,10]^{\top}$.

###  D.6 Ex211

The Ex211 problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-(42x_{1}-50y_{1}+44x_{2}+45x_{3}+47x_{4}+47.5x_{5}),$
$\displaystyle g_{1}(x,y)$ $\displaystyle=-(20x_{1}+y_{2}+4x_{5}-39).$
$\displaystyle
h(x)=\begin{bmatrix}x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}+x_{5}^{2}\\\
12x_{2}+11x_{3}+7x_{4}\end{bmatrix}.$ $\displaystyle 0\leq x_{i}\leq
1,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall
i=1,\ldots,5.$

The global maximum is $17$ with $x^{\star}=[1,1,0,1,0]^{\top}$.

###  D.7 Ex212

The Ex212 problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=10x_{6}+y_{1}+0.5(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}+x_{5}^{2}),$
$\displaystyle g_{1}(x,y)$
$\displaystyle=-(6x_{1}+3x_{2}+3x_{3}+2x_{4}+x_{5}-6.5),$ $\displaystyle
g_{2}(x,y)$ $\displaystyle=-(y_{2}-20).$ $\displaystyle
h(x)=\begin{bmatrix}10.5x_{1}+7.5x_{2}+3.5x_{3}+2.5x_{4}+1.5x_{5}\\\
10x_{1}+10x_{3}+x_{6}\end{bmatrix}.$ $\displaystyle 0\leq x_{i}\leq
30,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall
i=1,\ldots,6.$

The global maximum is $213$ with $x^{\star}=[0,1,0,1,1,20]^{\top}$.

###  D.8 g09

The g09 problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-y_{1}-x_{3}^{4}-3(x_{4}-11)^{2}-10x_{5}^{6}-7x_{6}^{2}-x_{7}^{4}+4x_{6}x_{7}+10x_{6}+8x_{7},$
$\displaystyle g_{1}(x,y)$ $\displaystyle=127-2x_{1}x_{2}-y_{2}-5x_{5},$
$\displaystyle g_{2}(x,y)$
$\displaystyle=282-7x_{1}-3x_{2}-10x_{3}^{2}-x_{4}+x_{5},$ $\displaystyle
g_{3}(x,y)$ $\displaystyle=196-23x_{1}+x_{2}^{2}-6x_{6}^{2}+8x_{7},$
$\displaystyle g_{4}(x,y)$
$\displaystyle=-4x_{1}^{2}-x_{2}^{2}+3x_{1}x_{2}-2x_{3}^{2}-5x_{6}+11x_{7}.$
$\displaystyle h(x)=\begin{bmatrix}(x_{1}-10)^{2}+5(x_{2}-12)^{2}\\\
3x_{2}^{4}+x_{3}+4x_{4}^{2}\end{bmatrix}.$ $\displaystyle-10\leq x_{i}\leq
10,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall
i=1,\ldots,7.$

The global maximum is $-680.63$ with
$x^{\star}=[2.33,1.95,-0.48,4.37,-0.62,1.04,1.59]^{\top}$.

###  D.9 Ex724

The Ex724 problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=-(y_{3}+0.4(x_{2}/x_{8})^{0.67}-x_{1}+10),$ $\displaystyle
g_{1}(x,y)$ $\displaystyle=-(0.0588x_{5}x_{7}+0.1x_{1}-1),$ $\displaystyle
g_{2}(x,y)$ $\displaystyle=-(0.0588x_{6}x_{8}+0.1x_{1}+0.1x_{2}-1),$
$\displaystyle g_{3}(x,y)$
$\displaystyle=-(4(x_{3}/x_{5})+2/y_{1}+0.0588(x_{7}/x_{3})^{1.3}-1),$
$\displaystyle g_{4}(x,y)$ $\displaystyle=-(y_{2}+0.0588x_{4}^{1.3}x_{8}-1).$
$\displaystyle h(x)=\begin{bmatrix}x_{3}^{0.71}x_{5}\\\
4(x_{4}/x_{6})+2/(x_{4}^{0.71}x_{6})\\\
0.4(x_{1}/x_{7})^{0.67}-x_{2}\end{bmatrix}.$ $\displaystyle 0.1\leq x_{i}\leq
10,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall
i=1,\ldots,8.$

The global maximum is $-3.92$ with
$x^{\star}=[6.35,2.34,0.67,0.53,5.95,5.32,1.04,0.42]^{\top}$.

###  D.10 Ex216

The Ex216 problem is defined by the following functions

$\displaystyle g_{0}(x,y)$
$\displaystyle=48x_{1}-0.5y_{1}-50x_{5}^{2}-50x_{6}^{2}-50x_{7}^{2}-50x_{8}^{2}-50x_{9}^{2}-50x_{10}^{2}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+42x_{2}+y_{3}+47x_{7}+42x_{8}+45x_{9}+46x_{10},$ $\displaystyle g_{1}(x,y)$
$\displaystyle=y_{2}-2x_{7}-6x_{8}-2x_{9}-2x_{1}0+4,$ $\displaystyle
g_{2}(x,y)$
$\displaystyle=6x_{1}-5x_{2}+8x_{3}-3x_{4}+x_{6}+3x_{7}+8x_{8}+9x_{9}-3x_{10}-22,$
$\displaystyle g_{3}(x,y)$
$\displaystyle=-5x_{1}+6x_{2}+5x_{3}+3x_{4}+8x_{5}-8x_{6}+9x_{7}+2x_{8}-9x_{10}+6,$
$\displaystyle g_{4}(x,y)$
$\displaystyle=y_{4}+3x_{7}-9x_{8}-9x_{9}-3x_{10}+23,$ $\displaystyle
g_{5}(x,y)$
$\displaystyle=-8x_{1}+7x_{2}-4x_{3}-5x_{4}-9x_{5}+x_{6}-7x_{7}-x_{8}+3x_{9}-2x_{10}+12.$
$\displaystyle
h(x)=\begin{bmatrix}100x_{1}^{2}+100x_{2}^{2}+100x_{3}^{2}+100x_{4}^{2}\\\
-2x_{1}6x_{2}-x_{3}-3x_{5}-3x_{6}\\\ 48x_{3}+45x_{4}+44x_{5}+41x_{6}\\\
9x_{1}+5x_{2}-9x_{4}+x_{5}-8x_{6}\end{bmatrix}.$ $\displaystyle 0\leq
x_{i}\leq 1,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\forall i=1,\ldots,10.$

The global maximum is $39$ with $x^{\star}=[1,0,0,1,1,1,0,1,1,1]^{\top}$.

## References

  * [1] D. Piga, M. Forgione, S. Formentin, A. Bemporad, Performance-oriented model learning for data-driven MPC design, IEEE control systems letters 3 (3) (2019) 577–582.
  * [2] J. A. Paulson, A. Mesbah, Data-driven scenario optimization for automated controller tuning with probabilistic performance guarantees, IEEE Control Systems Letters 5 (4) (2020) 1477–1482.
  * [3] F. Sorourifar, G. Makrygirgos, A. Mesbah, J. A. Paulson, A data-driven automatic tuning method for MPC under uncertainty using constrained Bayesian optimization, IFAC-PapersOnLine 54 (3) (2021) 243–250.
  * [4] E. A. del Rio Chanona, P. Petsagkourakis, E. Bradford, J. A. Graciano, B. Chachuat, Real-time optimization meets Bayesian optimization and derivative-free optimization: A tale of modifier adaptation, Computers & Chemical Engineering 147 (2021) 107249.
  * [5] J. Wu, S. Toscano-Palmerin, P. I. Frazier, A. G. Wilson, Practical multi-fidelity bayesian optimization for hyperparameter tuning, in: Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, Vol. 115 of Proceedings of Machine Learning Research, PMLR, 2020, pp. 788–798.
  * [6] I. Kapetanovic, Computer-aided drug discovery and development (CADDD): in silico-chemico-biological approach, Chemico-biological Interactions 171 (2) (2008) 165–176.
  * [7] S. Ju, T. Shiga, L. Feng, Z. Hou, K. Tsuda, J. Shiomi, Designing nanostructures for phonon transport via Bayesian optimization, Physical Review X 7 (2) (2017) 021024.
  * [8] A. M. Schweidtmann, A. D. Clayton, N. Holmes, E. Bradford, R. A. Bourne, A. A. Lapkin, Machine learning meets continuous flow chemistry: Automated optimization towards the Pareto front of multiple objectives, Chemical Engineering Journal 352 (2018) 277–282.
  * [9] J. Vrugt, J. W. Hopmans, J. Šimunek, Calibration of a two-dimensional root water uptake model, Soil Science Society of America Journal 65 (4) (2001) 1027–1037.
  * [10] L. Schultz, V. Sokolov, Bayesian optimization for transportation simulators, Procedia Computer Science 130 (2018) 973–978.
  * [11] J. Paulson, M. Martin-Casas, A. Mesbah, Fast uncertainty quantification in dynamic flux balance analysis models using sparse multi-element polynomial chaos, in: 2019 AIChE Annual Meeting, AIChE, 2019.
  * [12] L. M. Rios, N. V. Sahinidis, Derivative-free optimization: A review of algorithms and comparison of software implementations, Journal of Global Optimization 56 (3) (2013) 1247–1293.
  * [13] R. Eberhart, J. Kennedy, Particle swarm optimization, in: Proceedings of the IEEE International Conference on Neural Networks, Vol. 4, Citeseer, 1995, pp. 1942–1948.
  * [14] N. Hansen, S. D. Müller, P. Koumoutsakos, Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES), Evolutionary Computation 11 (1) (2003) 1–18.
  * [15] D. M. Mukhopadhyay, M. O. Balitanas, A. Farkhod, S.-H. Jeon, D. Bhattacharyya, Genetic algorithm: A tutorial review, International Journal of Grid and Distributed Computing 2 (3) (2009) 25–32.
  * [16] J. A. Nelder, R. Mead, A simplex method for function minimization, The Computer Journal 7 (4) (1965) 308–313.
  * [17] T. G. Kolda, R. M. Lewis, V. Torczon, Optimization by direct search: New perspectives on some classical and modern methods, SIAM Review 45 (3) (2003) 385–482.
  * [18] D. R. Jones, C. D. Perttunen, B. E. Stuckman, Lipschitzian optimization without the lipschitz constant, Journal of optimization Theory and Applications 79 (1) (1993) 157–181.
  * [19] R. Chen, M. Menickelly, K. Scheinberg, Stochastic optimization using a trust-region method and random models, Mathematical Programming 169 (2018) 447–487.
  * [20] F. E. Curtis, V. Kungurtsev, D. P. Robinson, Q. Wang, A stochastic-gradient-based interior-point algorithm for solving smooth bound-constrained optimization problems, arXiv preprint arXiv:2304.14907 (2023).
  * [21] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, N. De Freitas, Taking the human out of the loop: A review of Bayesian optimization, Proceedings of the IEEE 104 (1) (2015) 148–175.
  * [22] P. I. Frazier, A tutorial on Bayesian optimization, arXiv preprint arXiv:1807.02811 (2018).
  * [23] Y. Sui, A. Gotovos, J. Burdick, A. Krause, Safe exploration for optimization with Gaussian processes, in: International Conference on Machine Learning, PMLR, 2015, pp. 997–1005.
  * [24] D. Bergmann, K. Graichen, Safe Bayesian optimization under unknown constraints, in: In Proceedings of the Conference on Decision and Control, IEEE, 2020, pp. 3592–3597.
  * [25] F. Berkenkamp, A. Krause, A. P. Schoellig, Bayesian optimization with safety constraints: Safe and automatic parameter tuning in robotics, Machine Learning (2021) 1–35.
  * [26] D. Krishnamoorthy, F. J. Doyle, Safe Bayesian optimization using interior-point methods—applied to personalized insulin dose guidance, IEEE Control Systems Letters 6 (2022) 2834–2839.
  * [27] J. R. Gardner, M. J. Kusner, Z. E. Xu, K. Q. Weinberger, J. P. Cunningham, Bayesian optimization with inequality constraints., in: International Conference on Machine Learning, Vol. 2014, 2014, pp. 937–945.
  * [28] V. Picheny, R. B. Gramacy, S. Wild, S. L. Digabel, Bayesian optimization under mixed constraints with a slack-variable augmented Lagrangian, in: Proceedings of the 30th International Conference on Neural Information Processing Systems, 2016, pp. 1443–1451.
  * [29] M. J. Sasena, P. Papalambros, P. Goovaerts, Exploration of metamodeling sampling criteria for constrained global optimization, Engineering Optimization 34 (3) (2002) 263–278.
  * [30] R. Priem, N. Bartoli, Y. Diouane, A. Sgueglia, Upper trust bound feasibility criterion for mixed constrained bayesian optimization with application to aircraft design, Aerospace Science and Technology 105 (2020) 105980.
  * [31] C. Lu, J. A. Paulson, No-regret Bayesian optimization with unknown equality and inequality constraints using exact penalty functions, IFAC-PapersOnLine 55 (7) (2022) 895–902.
  * [32] W. Xu, Y. Jiang, C. N. Jones, Constrained efficient global optimization of expensive black-box functions, arXiv preprint arXiv:2211.00162 (2022).
  * [33] H. Chen, Lower rate of convergence for locating a maximum of a function, The Annals of Statistics (1988) 1330–1334.
  * [34] J. P. Eason, L. T. Biegler, A trust region filter method for glass box/black box optimization, AIChE Journal 62 (9) (2016) 3124–3136.
  * [35] B. Beykal, F. Boukouvala, C. A. Floudas, E. N. Pistikopoulos, Optimal design of energy systems using constrained grey-box multi-objective optimization, Computers & Chemical Engineering 116 (2018) 488–502.
  * [36] I. Bajaj, S. S. Iyer, M. F. Hasan, A trust region-based two phase algorithm for constrained black-box and grey-box optimization with infeasible initial point, Computers & Chemical Engineering 116 (2018) 306–321.
  * [37] S. H. Kim, F. Boukouvala, Surrogate-based optimization for mixed-integer nonlinear problems, Computers & Chemical Engineering 140 (2020) 106847.
  * [38] J. A. Paulson, C. Lu, COBALT: COnstrained Bayesian optimizAtion of computationaLly expensive grey-box models exploiting derivaTive information, Computers & Chemical Engineering 160 (2022) 107700.
  * [39] N. Srinivas, A. Krause, S. M. Kakade, M. Seeger, Gaussian process optimization in the bandit setting: No regret and experimental design, in: International Conference on Machine Learning, 2015, pp. 2171–2180.
  * [40] M. Blondel, O. Teboul, Q. Berthet, J. Djolonga, Fast differentiable sorting and ranking, in: International Conference on Machine Learning, PMLR, 2020, pp. 950–959.
  * [41] C. K. Williams, C. E. Rasmussen, Gaussian Processes for Machine Learning, Vol. 2, MIT Press Cambridge, MA, 2006.
  * [42] H. Liu, J. Cai, Y.-S. Ong, Remarks on multi-output Gaussian process regression, Knowledge-Based Systems 144 (2018) 102–121.
  * [43] D. Eriksson, M. Poloczek, Scalable constrained Bayesian optimization, in: International Conference on Artificial Intelligence and Statistics, PMLR, 2021, pp. 730–738.
  * [44] W. J. Maddox, M. Balandat, A. G. Wilson, E. Bakshy, Bayesian optimization with high-dimensional outputs, Advances in Neural Information Processing Systems 34 (2021) 19274–19287.
  * [45] R. Astudillo, P. I. Frazier, Thinking inside the box: A tutorial on grey-box Bayesian optimization, in: 2021 Winter Simulation Conference, IEEE, 2021, pp. 1–15.
  * [46] J. Močkus, On Bayesian methods for seeking the extremum, in: Optimization Techniques IFIP Technical Conference, Springer, 1975, pp. 400–404.
  * [47] D. R. Jones, M. Schonlau, W. J. Welch, Efficient global optimization of expensive black-box functions, Journal of Global Optimization 13 (4) (1998) 455\.
  * [48] J. T. Wilson, R. Moriconi, F. Hutter, M. P. Deisenroth, The reparameterization trick for acquisition functions, arXiv preprint arXiv:1712.00424 (2017).
  * [49] Y. Zhang, X. Zhang, P. Frazier, Constrained two-step look-ahead Bayesian optimization, Advances in Neural Information Processing Systems 34 (2021) 12563–12575.
  * [50] R. Astudillo, P. Frazier, Bayesian optimization of composite functions, in: International Conference on Machine Learning, PMLR, 2019, pp. 354–363.
  * [51] S. R. Chowdhury, A. Gopalan, On kernelized multi-armed bandits, in: International Conference on Machine Learning, PMLR, 2017, pp. 844–853.
  * [52] S. Vakili, K. Khezeli, V. Picheny, On information gain and regret bounds in Gaussian process bandits, in: International Conference on Artificial Intelligence and Statistics, PMLR, 2021, pp. 82–90.
  * [53] T. G. Epperly, E. N. Pistikopoulos, A reduced space branch and bound algorithm for global optimization, Journal of Global Optimization 11 (3) (1997) 287.
  * [54] V. J. Bowman, Permutation polyhedra, SIAM Journal on Applied Mathematics 22 (4) (1972) 580–589.
  * [55] K. Wang, B. Wilder, S.-c. Suen, B. Dilkina, M. Tambe, Improving GP-UCB algorithm by harnessing decomposed feedback, in: Machine Learning and Knowledge Discovery in Databases, Springer, 2020, pp. 555–569.
  * [56] M. Wüthrich, B. Schölkopf, A. Krause, Regret bounds for Gaussian-process optimization in large domains, Advances in Neural Information Processing Systems 34 (2021) 7385–7396.
  * [57] K. Chaloner, I. Verdinelli, Bayesian experimental design: A review, Statistical Science (1995) 273–304.
  * [58] M. Subotic, M. Tuba, N. Stanarevic, Different approaches in parallelization of the artificial bee colony algorithm, International Journal of Mathematical Models and Methods in Applied Sciences 5 (4) (2011) 755–762.
  * [59] M. Jamil, X.-S. Yang, A literature survey of benchmark functions for global optimisation problems, International Journal of Mathematical Modelling and Numerical Optimisation 4 (2) (2013) 150–194.
  * [60] L. A. Rastrigin, Systems of extremal control, Nauka (1974).
  * [61] S. Rahnamayan, H. R. Tizhoosh, M. M. Salama, A novel population initialization method for accelerating evolutionary algorithms, Computers & Mathematics with Applications 53 (10) (2007) 1605–1614.
  * [62] R. B. Gramacy, H. K. Lee, Cases for the nugget in modeling computer experiments, Statistics and Computing 22 (2012) 713–722.
  * [63] H. Rosenbrock, An automatic method for finding the greatest or least value of a function, The Computer Journal 3 (3) (1960) 175–184.
  * [64] R. Chelouah, P. Siarry, Tabu search applied to global optimization, European Journal of Operational Research 123 (2) (2000) 256–270.
  * [65] M. Laguna, R. Marti, Experimental testing of advanced scatter search designs for global optimization of multimodal functions, Journal of Global Optimization 33 (2005) 235–255.
  * [66] M. Styblinski, T.-S. Tang, Experiments in nonconvex optimization: stochastic approximation with function smoothing and simulated annealing, Neural Networks 3 (4) (1990) 467–483.
  * [67] M. S. Bazaraa, H. D. Sherali, C. M. Shetty, Nonlinear Programming: Theory and Algorithms, John wiley & sons, 2013.
  * [68] S. R. Nekoo, J. Á. Acosta, A. Ollero, A search algorithm for constrained engineering optimization and tuning the gains of controllers, Expert Systems with Applications 206 (2022) 117866.
  * [69] C. A. Floudas, P. M. Pardalos, C. Adjiman, W. R. Esposito, Z. H. Gümüs, S. T. Harding, J. L. Klepeis, C. A. Meyer, C. A. Schweiger, Handbook of test problems in local and global optimization, Vol. 33, Springer Science & Business Media, 2013.
  * [70] W. Hock, K. Schittkowski, Test examples for nonlinear programming codes, Journal of Optimization Theory and Applications 30 (1) (1980) 127–129.
  * [71] M. Tawarmalani, N. V. Sahinidis, Convexification and Global Optimization in Continuous and Mixed-integer Nonlinear Programming: Theory, Algorithms, Software, and Applications, Vol. 65, Springer Science & Business Media, 2013\.
  * [72] D. Karaboga, B. Akay, A modified artificial bee colony (ABC) algorithm for constrained optimization problems, Applied Soft Computing 11 (3) (2011) 3021–3031.
  * [73] W. Huyer, A. Neumaier, SNOBFIT–stable noisy optimization by branch and fit, ACM Transactions on Mathematical Software 35 (2) (2008) 1–25.
  * [74] scikit-quant (2023).
URL https://github.com/scikit-quant/scikit-quant

  * [75] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17 (2020) 261–272.  doi:10.1038/s41592-019-0686-2.
  * [76] M. J. Powell, The BOBYQA algorithm for bound constrained optimization without derivatives, Cambridge NA Report NA2009/06, University of Cambridge, Cambridge 26 (2009).
  * [77] C. Cartis, J. Fiala, B. Marteau, L. Roberts, Improving the flexibility and robustness of model-based derivative-free optimization solvers, ACM Transactions on Mathematical Software (TOMS) 45 (3) (2019) 1–41.
  * [78] N. Hansen, Y. Akimoto, P. Baudis, CMA-ES/pycma on Github, Zenodo, DOI:10.5281/zenodo.2559634 (Feb. 2019).  doi:10.5281/zenodo.2559634.
URL https://doi.org/10.5281/zenodo.2559634

  * [79] M. Balandat, B. Karrer, D. Jiang, S. Daulton, B. Letham, A. G. Wilson, E. Bakshy, BoTorch: A framework for efficient Monte-Carlo Bayesian optimization, Advances in Neural Information Processing Systems 33 (2020) 21524–21538.
  * [80] J. Larson, M. Menickelly, S. M. Wild, Derivative-free optimization methods, Acta Numerica 28 (2019) 287–404.
  * [81] C. Lu, J. A. Paulson, Constrained Upper Quantile Bound (2023).
URL https://github.com/PaulsonLab/Constrained_Upper_Quantile_Bound

  * [82] C. Zhu, R. H. Byrd, P. Lu, J. Nocedal, Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization, ACM Transactions on Mathematical Software 23 (4) (1997) 550–560.
  * [83] J. J. Moré, S. M. Wild, Benchmarking derivative-free optimization algorithms, SIAM Journal on Optimization 20 (1) (2009) 172–191.
  * [84] N. Bliznyuk, D. Ruppert, C. Shoemaker, R. Regis, S. Wild, P. Mugunthan, Bayesian calibration and uncertainty analysis for computationally expensive models using optimization and radial basis function approximation, Journal of Computational and Graphical Statistics 17 (2) (2008) 270–294.
  * [85] D. F. Mendoza, J. E. A. Graciano, F. dos Santos Liporace, G. A. C. Le Roux, Assessing the reliability of different real-time optimization methodologies, The Canadian Journal of Chemical Engineering 94 (3) (2016) 485–497.
  * [86] A. M. Schweidtmann, D. Bongartz, D. Grothe, T. Kerkenhoff, X. Lin, J. Najman, A. Mitsos, Deterministic global optimization with Gaussian processes embedded, Mathematical Programming Computation 13 (3) (2021) 553–581.