# The nature and composition of Jupiter’s building blocks derived from the
water abundance measurements by the Juno spacecraft

Olivier Mousis11affiliation: Aix Marseille Univ, CNRS, CNES, LAM, Marseille,
France<EMAIL_ADDRESS>, Jonathan I. Lunine22affiliation: Department of
Astronomy, Cornell University, Ithaca, NY 14853, USA , and Artyom
Aguichine11affiliation: Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
<EMAIL_ADDRESS>
###### Abstract

The microwave radiometer aboard the Juno spacecraft provided a measurement of
the water abundance found to range between $\sim$1 and 5.1 times the
protosolar abundance of oxygen in the near-equatorial region of Jupiter. Here,
we aim to combine this up-to-date oxygen determination, which is likely to be
more representative of the bulk abundance than the Galileo probe subsolar
value, with the other known measurements of elemental abundances in Jupiter,
to derive the formation conditions and initial composition of the building
blocks agglomerated by the growing planet, and that determine the heavy
element composition of its envelope. We investigate several cases of icy
solids formation in the protosolar nebula, from the condensation of pure ices
to the crystallization of mixtures of pure condensates and clathrates in
various proportions. Each of these cases correspond to a distinct solid
composition whose amount is adjusted in the envelope of Jupiter to match the O
abundance measured by Juno. The volatile enrichments can be matched by a wide
range of planetesimal compositions, from solids exclusively formed from pure
condensates or from nearly exclusively clathrates, the latter case providing a
slightly better fit. The total mass of volatiles needed in the envelope of
Jupiter to match the observed enrichments is within the $\sim$4.3–39
$M_{\oplus}$ range, depending on the crystallization scenario considered in
the protosolar nebula. A wide range of masses of heavy elements derived from
our fits is found compatible with the envelope’s metallicity calculated from
current interior models.

planets and satellites: individual (Jupiter) – planets and satellites: gaseous
planets – planets and satellites: interiors – planets and satellites:
formation – protoplanetary disks

##  1 Introduction

A long awaited measurement at Jupiter is the determination of its bulk oxygen
abundance. This measurement is key because it provides clues on the formation
conditions of the giant planet in the protosolar nebula (PSN) (Owen et al.,
1999; Gautier et al., 2001; Mousis et al., 2009, 2019). The in situ
measurements made by the Galileo probe up to 22 bar in Jupiter’s atmosphere
showed that the O abundance was subsolar ($\sim$0.46 times the protosolar O
abundance; Wong et al. (2004)). This result was mostly attributed to the
dynamics of the region within which the Galileo probe descended (Orton et al.,
1998), but alternative interpretations suggested that the measurement may
correspond to the bulk composition of the planet (Mousis et al., 2012).
However, the findings of the Galileo Probe were not conclusive because the
data showed that the abundance of water was still increasing with progressing
depth, until loss of signal. Recently, the microwave radiometer (MWR) aboard
the Juno spacecraft provided a measurement of the water abundance in the
equatorial region of Jupiter, from 0 to 4 degrees of north latitude, and in
the 0.7–30 bar pressure domain. The water abundance in this region of Jupiter
was found to range between $\sim$1 and 5 times the protosolar abundance of
oxygen, a value substantially higher than the Galileo probe determination, and
compatible with the presence of a moist adiabatic temperature profile at the
water condensation level (Li et al., 2020).

Here, we aim to combine this up-to-date oxygen determination, which is likely
to be more representative of the bulk abundance than the Galileo probe value,
with the other known measurements of elemental abundances in Jupiter, to
derive the formation conditions and initial composition of the building blocks
acquired by the growing planet, and that determined the composition of its
envelope, in the framework of the core accretion model (Pollack et al., 1996;
Inaba et al., 2003; Podolak et al., 2020). To do so, we also use the recent
Juno N abundance measurement that is found correlated with the O abundance in
Jupiter’s envelope (Li et al., 2020), as well as the in situ determinations of
Ar, Kr, Xe, C, N, and S by the Galileo probe (Mahaffy et al., 2000; Wong et
al., 2004), and the P abundance inferred from infrared observations by the
Cassini spacecraft (Fletcher et al., 2009) (see Table 1).

We investigate several cases of icy solids formation in the PSN, from the
condensation of pure ices to the crystallization of mixtures of pure
condensates and clathrates in various proportions. Each of these cases
corresponds to a distinct solid composition whose amount is adjusted in the
envelope of Jupiter to match the O abundance measured by Juno. The volatile
enrichments can be matched by a wide range of planetesimal compositions, from
solids exclusively formed from pure condensates or from nearly pure
clathrates, these latter providing the best fits. Estimates of the
corresponding amounts of heavy elements in the envelope are also provided and
compared to interior models derived from Juno observations.

Table 1: Elemental abundances and ratios to protosolar values in the upper troposphere of Jupiter Elements | Elemental abundances | Jupiter | References
---|---|---|---
|  | /Protosolar1 |
Ar/H | ($9.10\pm 1.80)\times 10^{-6}$ | $3.30\pm 0.65$ | Mahaffy et al. (2000)
Kr/H | ($4.65\pm 0.85)\times 10^{-9}$ | $2.38\pm 0.44$ | Mahaffy et al. (2000)
Xe/H | ($4.45\pm 0.85)\times 10^{-10}$ | $2.33\pm 0.45$ | Mahaffy et al. (2000)
C/H | ($1.19\pm 0.29)\times 10^{-3}$ | $4.01\pm 0.97$ | Wong et al. (2004)
N/H | ($2.04\pm 0.13)\times 10^{-4}$ | $2.75\pm 0.17$ | Bolton et al. (2017); Li et al. (2017, 2020)
O/H | ($2.45\pm 0.80)\times 10^{-4}$ | $0.46\pm 0.15$ | Wong et al. (2004)
| ($1.45^{+1.27}_{-0.93})\times 10^{-3}$ | $2.69^{+2.37}_{-1.72}$ | Li et al. (2020)
S/H | ($4.45\pm 1.05)\times 10^{-5}$ | $3.08\pm 0.73$ | Wong et al. (2004)
P/H | ($1.08\pm 0.06)\times 10^{-6}$ | $3.82\pm 0.20$ | Fletcher et al. (2009)

##  2 Building blocks composition model

In our model, the volatile phase incorporated in planetesimals is composed of
a mixture of pure ices, stoichiometric hydrates (such as NH3-H2O hydrate), and
clathrates that crystallized in the form of microscopic grains at various
temperatures in the outer part of the disk. Solids accreted by the growing
Jupiter are assumed to have preserved the volatile budget acquired during
condensation of the microscopic ices from which they assembled, regardless of
the way they evolved in the PSN (accretion of pebbles or large planetesimals,
products of collisions, etc).

We assume that the PSN is uniformly filled with H2O, CO, CO2, CH3OH, CH4, N2,
NH3, H2S, and PH3. These molecules are considered to be the dominant volatile
species in the PSN, assuming protosolar abundances for O, C, N, S, and P
(Asplund et al., 2009). Half of the sulfur is assumed to be in the form of H2S
with the other half forming refractory sulfide components (Pasek et al.,
2005), and all C forms CO, CO2, CH3OH and CH4, with the remaining O going into
H2O. We have set CO:CO2:CH3OH:CH4 = 10:30:1.67:1 in the gas phase of the disk.
The CO:CO2 ratio comes from ROSINA observations of the coma of comet
67P/Churyumov-Gerasimenko (hereafter 67P/C-G) (Läuter et al., 2019). The
CO:CH3OH:CH4 ratio is consistent with the production rates measured in the
southern hemisphere of 67P/C-G in October 2014 by the ROSINA instrument (Le
Roy et al., 2015). We also assume N2:NH3 = 1:1 in the nebula gas phase
(Fegley, 2000).

The process of volatile trapping in grains crystallized in the outer regions
of the PSN is determined by using the equilibrium curve of NH3 hydrate, those
of the various clathrates and pure condensates, as well as the thermodynamic
path detailing the evolution of temperature and pressure at 5.2 AU (i.e., the
current location of Jupiter) in the protoplanetary disk. We refer the reader
to the work of Mousis et al. (2009) for a full description of the calculation
method regarding the composition of the solids crystallized in the PSN. The
equilibrium curves of hydrates and clathrates derive from the compilation of
published experimental work by Lunine & Stevenson (1985), in which data are
available at relatively low temperatures and pressures. The equilibrium curves
of pure condensates used in our calculations derive from the compilation of
laboratory data from Lide (2002). The fits of the equilibrium data are taken
from Hersant et al. (2004) and Mousis et al. (2008).

The disk model employed is the one described in Aguichine et al. (2020) and
Mousis et al. (2020), to which the reader is referred for details. In a few
words, our time-dependent PSN model is governed by the following differential
equation (Lynden-Bell & Pringle, 1974):

$\displaystyle\frac{\partial\Sigma_{\mathrm{g}}}{\partial
t}=\frac{3}{r}\frac{\partial}{\partial r}\left[r^{1/2}\frac{\partial}{\partial
r}\left(r^{1/2}\Sigma_{\mathrm{g}}\nu\right)\right].$ (1)

This equation describes the time evolution of a viscous accretion disk of
surface density $\Sigma_{\mathrm{g}}$ of viscosity $\nu$, assuming hydrostatic
equilibrium in the $z$ direction. The viscosity $\nu$ is calculated in the
framework of the $\alpha$-formalism (Shakura & Sunyaev, 1973). In our
approach, the midplane temperature $T$ of the disk is expressed as the sum of
viscous heating and star background irradiation $T_{\mathrm{amb}}$= 10 K:

$T^{4}=\frac{1}{2\sigma_{\mathrm{SB}}}\left(\frac{3}{8}\tau_{\mathrm{R}}+\frac{1}{2\tau_{\mathrm{P}}}\right)\Sigma_{\mathrm{g}}\nu\Omega_{\mathrm{K}}^{2}+T^{4}_{\mathrm{amb}},$
(2)

where $\sigma_{\mathrm{sb}}$ is the Stefan-Boltzmann constant,
$\Omega_{\mathrm{K}}$ is the keplerian frequency, $\tau_{\mathrm{R}}$ and
$\tau_{\mathrm{P}}$ are the Rosseland and Planck optical depth, respectively
(see Aguichine et al. (2020) and Mousis et al. (2020) for details). Only
viscous heating is considered as an energy source in our model, allowing the
disk temperature to decrease down to $\sim$10 K in the outer PSN. Here viscous
means eddy viscosity, namely that provided by eddies which transport mass and
heat by bulk motions, since molecular viscosity is very small (Stevenson,
1990). This corresponds to the case where the outer parts of the disk are
protected from solar irradiation by shadowing effect of the inner disk parts.
The evolution of the disk starts with an initial profile given by
$\Sigma_{\mathrm{g}}\nu\propto\exp\left(-r^{2-p}\right)$, with $p=\frac{3}{2}$
for an early disk (Lynden-Bell & Pringle, 1974). In our computations, the
initial disk mass is fixed to 0.1 M⊙. The computational box is set equal to
500 AU, allowing 99% of the disk mass to be encapsulated within $\sim$100 AU.
The disk initial mass accretion rate onto the Sun is set to $10^{-7.6}$ M⊙
yr-1 (Hartmann et al., 1998), and the viscosity parameter $\alpha$ is fixed
equal to 3 $\times$ 10-3, which is well within the 10-4–10-2 range commonly
adopted for PSN models (Mousis et al., 2020).

In the following, we assume that condensate or clathrate grains quickly
decouple from gas, due to growth and planetesimal formation. The enrichment of
species $X$ with respect to its protosolar abundance is then given by:

$\frac{X/H_{2}(t)}{X/H_{2}|_{\odot}}=\frac{\Sigma_{\mathrm{g}}(R,t_{X})}{\Sigma_{\mathrm{g}}(R,t)},$
(3)

where if $t_{X}$ is the time at which species $X$ is trapped or condensed at
given distance $r$ in the PSN, This relation is valid for the determination of
the volatile enrichments in Jupiter’s feeding zone, and subsequently in its
envelope. A set of volatile enrichments is simultaneously calculated in
Jupiter when Eq. 3 is fitted to the abundance measurement of a given element
in the envelope.

Figure 1 represents the condensation sequence of the different volatiles at
5.2 AU in the PSN, for various initial H2O abundances. The domain of stability
of each ice considered is the region located below its corresponding
equilibrium curve. In all panels, it is assumed that CO:CO2:CH3OH:CH4 =
10:30:1.67:1, H2S/H2 = 0.5 $\times$ (S/H2)☉, and N2:NH3 = 1:1 in the PSN gas
phase, based on the aforementioned assumptions. Four cases are investigated.
Case 1 (panel (a)) corresponds to the case where the abundances of all
elements are protosolar (Asplund et al., 2009), assuming volatile species form
only pure condensates in the PSN. Case 2 (panel (b)) corresponds to the same
assumptions as in case 1 regarding the PSN gas phase composition, but here it
is assumed that all available water is used to form hydrates in the PSN. In
this situation, NH3-H2O hydrate forms first with decreasing disk temperature,
and then almost 75% of available H2S is trapped in clathrates. The remaining
species subsequently form pure condensates at lower disk temperatures. Note
that CO2 and CH3OH are both represented in their pure ice forms in all our
condensation sequences. Pure CO2 is the only condensate that crystallizes at
higher temperature than its associated clathrate in this pressure range and no
experimental data concerning the equilibrium curve of CH3OH clathrate has been
reported in the literature. Case 3 (panel (c)) corresponds to the same
assumptions as in case 2, except that O/H2 is 1.5 $\times$ (O/H2)☉. Here,
NH3-H2O hydrate forms, and H2S, PH3, Xe, CH4, and CO are fully trapped in
clathrates. No more than $\sim$67% of available N2 in enclathrated because of
the lack of available free water, and the remaining volatile species form pure
condensates at lower disk temperatures. Case 4 (panel (d)) corresponds to the
same assumptions as in case 3, except that O/H2 is 2 $\times$ (O/H2)☉. Here,
all volatiles except CH3OH and CO2, form hydrates in the PSN. Note that we do
not consider mixed clathrates, that is, incorporation of species at
temperatures other than where their pure clathrates form. This does not
significantly affect the final result in terms of comparison with Jupiter
data.

Figure 2 shows the composition of the volatile phase embedded in solids formed
in the PSN resulting from our model. Calculations represented in panels (a–d)
follow the cases 1–4 depicted above, respectively. CO2 is by far the main
volatile species in panels (a) and (b) while H2O dominates in panels (c) and
(d). The high abundance of H2O in those panels results from the choice of a
larger O abundance in the PSN. In panels (b–d), the icy part of planetary
building blocks is essentially made of a mixture of pure condensates and
clathrates in proportions fixed by the PSN oxygen abundance.

Figure 1:  Formation conditions of icy solids in the PSN defined by the
intersection between various equilibrium curves and the disk cooling curve at
5.2 AU. Arrows pointing down in the cooling curves indicate the direction of
time evolution. Species remain in the gas phase above the equilibrium curves.
Below, they are trapped as clathrates or simply condense. Panel (a):
equilibrium curves of pure condensates (red lines). Abundances of various
elements are protosolar (Asplund et al., 2009), with molecular mixing ratios
given in the text. Panel (b): same gas phase conditions as in panel (a) but
with the assumption that all available H2O is used to form NH3 hydrate
(NH3-H2O) and H2S clathrate (H2S-5.75H2O) (blue lines) in the PSN. Remaining
species form pure condensates (red lines). Panel (c): same as panel (b),
except that O/H2 is 1.5 $\times$ (O/H2)☉ in the PSN, implying that the number
of entrapped species is more important. Panel (d): same as panels (b) and (c),
except that O/H2 is 2 $\times$ (O/H2)☉ in the PSN, implying that all species
potentially forming clathrates become entrapped.

Figure 2: Composition of the volatile phase incorporated in solids in the
cases the volatile part of the building blocks is formed from pure condensates
only (panel (a)), a mixture of pure condensates, NH3 hydrate and various
clathrates (panels (b) and (c)), and from NH3 hydrate and clathrates only
(panel (d)). Initial gas phase conditions in panels (a)–(d) are those depicted
in the text and given in Fig. 1.

##  3 Results

Figure 3 represents the fits of the volatile enrichments (relative to
protosolar abundances) observed in Jupiter, and summarized in Table 1. Panels
(a–d) are calculated following the conditions depicted in cases 1–4,
respectively. All panels have been fitted to the oxygen abundance measured by
Juno in Jupiter’s atmosphere. The fits consist in reproducing the mass range
of oxygen measured in Jupiter by adjusting the amount of O–bearing solids
(here H2O, CO, CO2, and CH3OH) formed in the PSN and that is injected in the
envelope, in all cases of planetesimals compositions. All our calculations
assume the envelope is homogeneously mixed. Panels (a–d) show that C, N, P,
Kr, and Xe elemental enrichments are always fitted by our different
compositions models of the building blocks accreted by Jupiter. S is fitted in
panels (a–b), which both correspond to the assumption that all elements,
including O, are in protosolar abundances in the initial gas phase of the PSN.
Panels (c–d), which assume both supersolar O abundances and clathrate
formation in Jupiter’s feeding zone, present calculations of S enrichments
slightly lower than the measured value. However, an increase by $\sim$6–12$\%$
of the assumed H2S abundance in the PSN, this latter being poorly assessed
(Pasek et al., 2005), or the delivery to Jupiter of other S–bearing solids
such as SO2 in the volatile phase or FeS in the mineral phase, would easily
enable the matching between the measured and calculated enrichments. The Ar
enrichment is only marginally fitted in panel (d), which corresponds to case 4
(full clathration and O/H2 = 2 $\times$ (O/H2)☉ in the PSN). Only the full
trapping of Ar in clathrates at higher temperature than its condensation in
the PSN (see Fig. 1) allows the matching of the observed enrichment.

The fitted enrichments also translate into a mass range of icy solids
incorporated in Jupiter’s envelope following the compositions shown in Fig. 2.
Table 2 summarizes these mass ranges, assuming an average envelope mass of
$\sim$301.8 $M_{\oplus}$ (Wahl et al., 2017). The mass ranges calculated for
cases 1–4 must be seen as minimums because the accreted solids should harbor a
significant fraction of minerals and metals. Assuming a fraction of rocks
(minerals + metals) half the total mass of condensed matter in the PSN for a
gas of protosolar composition (Lodders, 2003), the mass of accreted solids
should be then twice the values indicated for cases 1 and 2 in Table 2, giving
a total mass range between 8.6 and 45.2 $M_{\oplus}$. The same extra mass of
rocks, i.e. between $\sim$4.3 and 22.6 $M_{\oplus}$ should be added to cases 3
and 4 since the extra oxygen added to the PSN gas in our model is assumed to
go only with the volatiles.

Table 2: Mass of volatiles needed in the envelope of Jupiter to match the observed enrichments | Mass of H2O ice ($M_{\oplus}$) | Mass of ices ($M_{\oplus}$)
---|---|---
Case 1 | 0.4–1.9 | 4.3–22.6
Case 2 | 0.4–1.9 | 4.3–22.7
Case 3 | 2.3–12.2 | 5.9–30.7
Case 4 | 4.1–21.5 | 7.4–39.0

Figure 3: Ratio of Jovian to protosolar abundances in the cases the volatile
part of the building blocks is formed from pure condensates only (panel (a)),
a mixture of pure condensates, NH3 hydrate and various clathrates (panels (b)
and (c)), and from NH3 hydrate and clathrates only (panel (d)). Initial gas
phase conditions in panels (a)–(d) are those given in Fig. 1. Blue and red
bars correspond to observations and model, respectively. Initial gas phase
conditions in panels (a)–(d) are those depicted in the text and given in Fig.
1. The green arrow pointing up illustrates that fact that the calculated S
enrichment is potentially higher if one assumes H2S/H2 $>$ 0.5 $\times$
(S/H2)☉ in the PSN gas phase.

##  4 Conclusion

Our work shows that a wide range of planetesimal compositions allows the fits
of the volatiles enrichments observed in Jupiter’s envelope, including the
recent Juno O determination. Ices either formed from pure condensates or from
clathrates can match the observed enrichments in a satisfying manner. We find
however that the fits of the volatile enrichments in the case of planetesimals
essentially formed from clathrates is slightly more compelling since it allows
marginally the matching of the Ar abundance in Jupiter, contrary to the case
of pure condensates. The agglomeration of planetesimals from clathrates
classically requires a supersolar abundance of oxygen at their formation
location (Gautier et al., 2001; Mousis et al., 2009). Such an effect in the
PSN can occur at the location of the snowline, where the outward diffusion of
vapor increases the local abundance of solid water (Stevenson & Lunine, 1988;
Cyr et al., 1999; Mousis et al., 2019).

The mass range of volatiles needed to be injected in solid form in Jupiter’s
envelope is within $\sim$4.3–39 $M_{\oplus}$, if one considers all
investigated cases. These values are well within the mass range of heavy
elements derived from Juno’s preliminary measurements of gravitational moments
(Wahl et al., 2017). Based on the interior model of Militzer & Hubbard (2013),
Wahl et al. (2017) showed that the mass of heavy elements is within the
$\sim$3.3–15.9 $M_{\oplus}$ range, depending on the type of core, which can be
diluted or compact. Using the model of Becker et al. (2013), these authors
also derive a larger amount of heavy elements in the envelope, which is within
the $\sim$14.5–40 $M_{\oplus}$ range. Even if a substantial fraction of
refractory component is taken into account in our calculations of planetesimal
compositions, a wide range of masses of heavy elements derived from our fits
is compatible with the envelope’s metallicity calculated from current interior
models.

Our calculations of the volatile budget in the envelope of Jupiter assume the
envelope is homogeneously mixed. This may be inaccurate, as some recent models
suggest a composition gradient in the envelope (Leconte & Chabrier, 2012;
Debras & Chabrier, 2019). However, further theoretical and experimental
studies of hydrogen-helium mixtures are needed below $\sim$100 GPa to assess
these hypotheses, a region where the models are critically sensitive to
changes in the equation of state (Wahl et al., 2017). Also, one cannot exclude
that the observed volatile enrichments correspond to a late planetesimal
accretion which, coupled with a rapid infall of gas, could have led to a
mixing of accreted material throughout the outer regions, which may explain
the supersolar metallicity (Podolak et al., 2020). However, in this scenario,
our results regarding the composition and structure of the volatile phase
embedded in planetesimals would remain valid.

Note that the Juno O determination in Jupiter considered in our work is
derived from a 1–$\sigma$ uncertainty. A 2–$\sigma$ uncertainty allows a H2O
abundance to be as low as 0.1 times the solar abundance (Li et al., 2020), a
value lower than the one ($\sim$0.5 times solar) inferred by the Galileo probe
(Wong et al., 2004). If the bulk O abundance in Jupiter is significantly
subsolar, this implies that the planet could belong to the category of carbon-
rich planets, i.e. those with C/O $\geq$ 1 in their envelopes (Mousis et al.,
2012). This would imply that water ice was heterogeneously distributed over
several AU beyond the snow line in the protosolar nebula and that the fraction
of water contained in icy planetesimals was a strong function of their
formation location and time.

Interestingly, for the sake of comparison, we have also used the compilation
of protosolar abundances derived from Lodders et al. (2009). Despite the
important deviations noted for the C/O and Ar/O ratios in Lodders et al.
(2009), which are +15% and -17% compared with Asplund et al. (2009),
respectively, the fits of the Jovian abundances with this database lead to
similar conclusions.

One should recall that the present study is based on the Juno MWR measurements
performed only at the equator and that did not penetrate deeper than the
$\sim$30 bar region, without using the data from the longest wavelength, 0.6
GHz channel. Future Juno measurements at other latitudes, combined with
updated interior structure models, should decrease the size of the O error bar
and enable our model to provide more stringent insights on Jupiter’s formation
conditions in the PSN. Finally, our calculations assume that the chemical
composition of the disk is homogeneous in the formation region of Jupiter’s
building blocks beyond the snowline. This assumption does not consider the
variations in the abundances of the various ices at the locations of their
condensation lines (Mousis et al., 2020, 2021). Future models describing the
composition of the solids formed in the PSN will have to consider these
physical processes.

O.M. acknowledges support from CNES. J.L. was supported by the Juno project
through a contract with Southwest Research Institute. We thank the anonymous
referee for the helpful comments, and constructive remarks on this manuscript.

## References

  * Aguichine et al. (2020) Aguichine, A., Mousis, O., Devouard, B., et al. 2020, ApJ, 901, 97. doi:10.3847/1538-4357/abaf47
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., et al. 2009, ARA&A, 47, 481. doi:10.1146/annurev.astro.46.060407.145222
  * Becker et al. (2013) Becker, A., Nettelmann, N., Holst, B., et al. 2013, Phys. Rev. B, 88, 045122. doi:10.1103/PhysRevB.88.045122
  * Bolton et al. (2017) Bolton, S. J., Adriani, A., Adumitroaie, V., et al. 2017, Science, 356, 821. doi:10.1126/science.aal2108
  * Cyr et al. (1999) Cyr, K. E., Sharp, C. M., & Lunine, J. I. 1999, J. Geophys. Res., 104, 19003. doi:10.1029/1998JE000617
  * Debras & Chabrier (2019) Debras, F. & Chabrier, G. 2019, ApJ, 872, 100. doi:10.3847/1538-4357/aaff65
  * Fegley (2000) Fegley, B. 2000, Space Sci. Rev., 92, 177. doi:10.1023/A:1005286910756
  * Fletcher et al. (2009) Fletcher, L. N., Orton, G. S., Teanby, N. A., et al. 2009, Icarus, 202, 543. doi:10.1016/j.icarus.2009.03.023
  * Gautier et al. (2001) Gautier, D., Hersant, F., Mousis, O., et al. 2001, ApJ, 550, L227. doi:10.1086/319648
  * Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., et al. 1998, ApJ, 495, 385
  * Hersant et al. (2004) Hersant, F., Gautier, D., & Lunine, J. I. 2004, Planet. Space Sci., 52, 623. doi:10.1016/j.pss.2003.12.011
  * Inaba et al. (2003) Inaba, S., Wetherill, G. W., & Ikoma, M. 2003, Icarus, 166, 46. doi:10.1016/j.icarus.2003.08.001
  * Läuter et al. (2019) Läuter, M., Kramer, T., Rubin, M., et al. 2019, MNRAS, 483, 852. doi:10.1093/mnras/sty3103
  * Leconte & Chabrier (2012) Leconte, J. & Chabrier, G. 2012, A&A, 540, A20. doi:10.1051/0004-6361/201117595
  * Le Roy et al. (2015) Le Roy, L., Altwegg, K., Balsiger, H., et al. 2015, A&A, 583, A1. doi:10.1051/0004-6361/201526450
  * Li et al. (2020) Li, C., Ingersoll, A., Bolton, S., et al. 2020, Nature Astronomy, 4, 609. doi:10.1038/s41550-020-1009-3
  * Li et al. (2017) Li, C., Ingersoll, A., Janssen, M., et al. 2017, Geophys. Res. Lett., 44, 5317. doi:10.1002/2017GL073159
  * Lide (2002) Lide, D. R. 2002, CRC Handbook of chemistry and physics : a ready-reference book of chemical and physical data, 83rd ed., by David R. Lide. Boca Raton: CRC Press, ISBN 0849304830, 2002
  * Lodders et al. (2009) Lodders, K., Palme, H., & Gail, H.-P. 2009, Landolt Börnstein, 4B, 712. doi:10.1007/978-3-540-88055-4$\\_$34
  * Lodders (2003) Lodders, K. 2003, ApJ, 591, 1220. doi:10.1086/375492
  * Lunine & Stevenson (1985) Lunine, J. I. & Stevenson, D. J. 1985, ApJS, 58, 493. doi:10.1086/191050
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603
  * Mahaffy et al. (2000) Mahaffy, P. R., Niemann, H. B., Alpert, A., et al. 2000, J. Geophys. Res., 105, 15061. doi:10.1029/1999JE001224
  * Militzer & Hubbard (2013) Militzer, B. & Hubbard, W. B. 2013, ApJ, 774, 148. doi:10.1088/0004-637X/774/2/148
  * Mousis et al. (2021) Mousis, O., Aguichine, A., Bouquet, A., et al. 2021, PSJ, 2, 72. doi:10.3847/PSJ/abeaa7
  * Mousis et al. (2020) Mousis, O., Aguichine, A., Helled, R., et al. 2020, Phil. Trans. R. Soc. A, 378: 20200107. http://dx.doi.org/10.1098/rsta.2020.0107
  * Mousis et al. (2019) Mousis, O., Ronnet, T., & Lunine, J. I. 2019, ApJ, 875, 9. doi:10.3847/1538-4357/ab0a72
  * Mousis et al. (2012) Mousis, O., Lunine, J. I., Madhusudhan, N., et al. 2012, ApJ, 751, L7. doi:10.1088/2041-8205/751/1/L7
  * Mousis et al. (2009) Mousis, O., Marboeuf, U., Lunine, J. I., et al. 2009, ApJ, 696, 1348. doi:10.1088/0004-637X/696/2/1348
  * Mousis et al. (2008) Mousis, O., Alibert, Y., Hestroffer, D., et al. 2008, MNRAS, 383, 1269. doi:10.1111/j.1365-2966.2007.12653.x
  * Orton et al. (1998) Orton, G. S., Fisher, B. M., Baines, K. H., et al. 1998, J. Geophys. Res., 103, 22791. doi:10.1029/98JE02380
  * Owen et al. (1999) Owen, T., Mahaffy, P., Niemann, H. B., et al. 1999, Nature, 402, 269. doi:10.1038/46232
  * Pasek et al. (2005) Pasek, M. A., Milsom, J. A., Ciesla, F. J., et al. 2005, Icarus, 175, 1. doi:10.1016/j.icarus.2004.10.012
  * Podolak et al. (2020) Podolak, M., Haghighipour, N., Bodenheimer, P., et al. 2020, ApJ, 899, 45. doi:10.3847/1538-4357/ab9ec1
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62. doi:10.1006/icar.1996.0190
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 500, 33
  * Stevenson (1990) Stevenson, D. J. 1990, ApJ, 348, 730. doi:10.1086/168282
  * Stevenson & Lunine (1988) Stevenson, D. J. & Lunine, J. I. 1988, Icarus, 75, 146. doi:10.1016/0019-1035(88)90133-9
  * Wahl et al. (2017) Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649. doi:10.1002/2017GL073160
  * Wong et al. (2004) Wong, M. H., Mahaffy, P. R., Atreya, S. K., et al. 2004, Icarus, 171, 153. doi:10.1016/j.icarus.2004.04.010


Send correspondence to<EMAIL_ADDRESS>
# Towards new servo control algorithms at the TNG telescope

P. Schipani  INAF - Osservatorio Astronomico di Capodimonte, Salita Moiariello
16, Naples, Italy  M. Gonzalez  INAF - FGG, TNG, Rambla J.A. Fernández Pérez
7, Breña Baja (TF), Spain  F. Perrotta  INAF - Osservatorio Astronomico di
Capodimonte, Salita Moiariello 16, Naples, Italy  S. Savarese  INAF -
Osservatorio Astronomico di Capodimonte, Salita Moiariello 16, Naples, Italy
M. Colapietro  INAF - Osservatorio Astronomico di Capodimonte, Salita
Moiariello 16, Naples, Italy  A. Ghedina  INAF - FGG, TNG, Rambla J.A.
Fernández Pérez 7, Breña Baja (TF), Spain  M. Hernandez Diaz  INAF - FGG, TNG,
Rambla J.A. Fernández Pérez 7, Breña Baja (TF), Spain  H. Ventura  INAF - FGG,
TNG, Rambla J.A. Fernández Pérez 7, Breña Baja (TF), Spain

###### Abstract

The servo control algorithms of the TNG, developed in the nineties, have been
working for more than 20 years with no major updates. The original hardware
was based on a VME-bus based platform running a real time operating system, a
rather popular choice for similar applications at the time. Recently, the
obsolescence of the hardware and the lack of spares pushed the observatory
towards a complete replacement of the electronics that is now being
implemented in steps, respecting the basic requirement of never stopping the
observatory night operations. Within the framework of this major hardware
work, we are taking the opportunity to review and update the existing control
schemes. This servo control update, crucial for the telescope performance,
envisages a new study from scratch of the controlled plant, including a re-
identification of the main axes transfer functions and a re-design of the
control filters in the two nested position and speed loops. The ongoing work
is described, including preliminary results in the case study of the azimuth
axis and our plans for possible further improvements.

###### keywords:

TNG, Telescope, Tracking, Servo Control

##  1 INTRODUCTION

The 3.5-m TNG telescope saw first light in La Palma (Canary Islands) at Roque
de Los Muchachos observatory in 1998. Its axes control system had been
realized few years before: all the control hardware had been installed in
1995, completing the procurement of components and custom parts during the
preliminary integration in Italy. In the original design, the speed loop
control had been implemented adopting a custom hardware developed by the
Sierracin Magnedyne company, including motors, power amplifiers and control
CPUs, whereas the position loop had been realized in-house by the Italian
community, using all commercial components based on the VME bus technology and
a real-time operating system: this system has been working every night for
more than 20 years.

Only few relatively minor interventions have been necessary over the years,
but after more than two decades, the observatory was running out of spare
parts. Soon after the telescope realization, the stock could no longer be
updated for the custom parts. Of course, over the years this became a critical
issue because of aging and obsolescence, including the commercial components
as well. Thus, recently the observatory started a project to replace all the
control hardware components but motors and feedback sensors: this work is
described in Ref. 1. Obviously, this kind of modernization efforts are common
for telescopes operating since a few decades (e.g. see Ref. 2, 3, 4).

Within this effort, it was decided to maintain the same software, which had
proved to be effective and reliable for many years, so avoiding a full
recommissioning of the system. Thus, the control software has been ported to
the new platform with a heavy and accurate work of translation from the old
real-time textual environment to the National Instruments LabVIEW graphical
programming language. Thus, the telescope is still running the “same” axes
control software albeit under a totally different environment.

As the whole intervention involves critical system components (e.g. the axes
drives), possible effects on the axes transfer function, and consequently on
the servo control performance, can be expected. This motivated the team to
plan a complete check and redesign of the servo control algorithms, under the
strict requirement of never stopping the telescope operations. This is the
main objective of this paper, which describes the case of the azimuth axis.

---
Figure 1:  Open loop model identification. The white noise is injected into
the torque input, the differentiated encoder signal is recorded synchronously
as system output. The I/O sample pairs are analyzed in Matlab through a laptop
computer.
---
Figure 2:  Open-loop transfer functions of the azimuth axis at three different
altitude angles: ALT=90∘ (Zenith), ALT=50∘, ALT=20∘ (observational limit).
---
Figure 3:  Open-loop transfer functions of the azimuth axis at ALT=50∘. Blue:
raw noisy data; Red: 50-state identified model; Yellow: 14-state model after
Hankel reduction.

##  2 System identification

The first step of the work has been the identification of the azimuth model.
The derivation of the transfer function has been done purely from field test
data, as no analytical model was available for comparison. A white noise
signal has been applied at the torque input to excite the dynamics of the
telescope. The white noise is a random signal of even spectrum, able to excite
evenly each frequency component of the telescope dynamics. The corresponding
velocity output has been measured by differentiating the position encoder
output.

It is worth noticing that, in the original implementation, the velocity
feedback was provided by custom tachometers; after several years they started
to fail and were replaced by an encoder based system, providing a back-
compatible analogue speed signal. Finally, in the new system it has been
decided to work with “software tachometers” obtained by differentiating the
encoder readings. This is the same solution adopted, e.g. at the VLT (Ref. 5)
and the VST (Ref. 6).

To be effective, the white noise procedure has been tuned by varying
parameters like the amplitude of the input noise and the duration of the test.
After some tuning, the test was run at 500Hz frequency for an overall duration
of 5 minutes, corresponding to 150,000 samples of input-output pairs. Longer
acquisition times did not produce any difference. The amplitude of the input
signal was determined after a trade-off between the opposite needs of exciting
sufficiently the axis dynamics and prevent damages to the system. The average
input signal was null, so the axis remained always approximately in the same
position.

The test configuration was extremely easy, as no external equipment was
adopted. Everything was based on dedicated software features prepared in
advance, fully integrated into the new graphical implementation of the control
software, adopting a post-processing code developed in Matlab. The white noise
input was purely digital as well as the velocity output. The sample vectors
were analyzed in quasi-real-time within the Matlab environment. Figure 1 shows
the basic diagram of the open-loop system identification.

In the case of azimuth, it is well known that the transfer function depends on
the altitude position (Ref. 7). The model identification has been repeated at
the two extremes of the altitude range (ALT=20∘ and ALT=90∘) and in an
intermediate position (ALT=50∘). As expected, the three transfer functions are
similar, but not identical (Figure 2). The result at ALT=50∘ is assumed to be
a good compromise, representing a first order approximation of the axis model.

In principle, the identified model might be compared with an analytical model
coming from FEA data, which would help to identify the physical meaning of all
resonances in the model. Nevertheless, an accurate FEA model is not available
for the TNG and any old simplified model, dating back to design time, could
not be deemed reliable. Also, the purpose of this work is the design of the
servo control for an existing telescope, rather than the study of a new
optomechanical system still in design phase. Consequently, all the work was
based uniquely on the identified transfer function. The identified model has
been the starting point for the servo control design, described in the
following sections.

It is worth mentioning that for some control strategies (e.g. model-based
control schemes) an accurate identified model is mandatory (Ref. 7). So far,
the old servo control architecture has been kept for this work, but a model-
based controller would be a viable alternative. Working on a model-based
controller (e.g. LQG) design, it is generally necessary to reduce the order of
the model to get an accurate, but simple, representation of the system. The
order of the system must be decreased saving only an appropriate number of
modes that still represent the system with good approximation. All the others,
adding unuseful complexity to the controller design, shall be excluded.

In the case under study, the azimuth axis was identified with a model of order
50, i.e. the dynamical matrix of the space-state representation is 50x50. This
order is too high for the implementation of e.g. a LQG controller and likely
includes dynamics that do not represent the system, possibly measurement
noise. Different reduction techniques exist (Ref. 8); here the Hankel norm
technique has been adopted, computing the Hankel singular vaues of the
dynamical system and removing the modes with the smallest norms, under a given
threshold. The threshold must be tuned in order to keep an adequate
representation of the system. Figure 3 compares raw data, the 50-state azimuth
model and a balanced 14-state model reduced by truncation of the smallest
Hankel norms; here the 50- and 14-state model differs significantly only for
the high frequency dynamics, above 100Hz.

---
Figure 4:  The speed loop scheme. The inner structure of the speed controller
is shown as a block diagram. On top, the Bode diagrams of the filters.
---
Figure 5:  Example of a velocity error time series, for a speed command of 80
arcsec/s (RMS=0.4 arcsec/s).

##  3 Speed Loop Control

The telescope main axes adopt a two-loop control structure with a speed loop
nested inside the position loop. Adopting this structure, the problem is
naturally decoupled in two steps, to be performed in the right order:

  * •

design of the speed controller, which defines the closed-loop speed transfer
function

  * •

design of the position controller

Thus, the speed controller design is the first step. A good speed control is
necessary, but not sufficient, for a good star tracking. When the velocity
control works sufficiently well, the control engineer can start dealing with
the position control loop.

The implementation of a series of second order filters offers a full
flexibility for the speed controller design. They have an immediate
interpretation, e.g. they can be low-pass, Notch, PID filters. The serial
connection of a number of simple blocks allows to shape any possible transfer
function in the frequency domain. Here, we preliminarily estimated that a
series of four filters could give enough degrees of freedom. However, a good
solution was identified using just three of them. The filters have been
designed in the frequency domain, by placing poles and zeros where
appropriate. In the current implementation, the filter no.1 and no.2 are
identical low-pass filters and the filter no.3 is a PI controller. An external
gain has been modulated in order to increase the bandwidth of the closed speed
loop as much as possible, preserving simultaneously the system stability.
Figure 4 shows the speed closed-loop with a schematic description of the speed
controller. Figure 5 reports an example of speed error time series during the
tuning phase.

##  4 Position Loop Control

The fundamental goal for a telescope axis servo control is the quality of
tracking. The control loops must be designed carefully to constrain the
tracking error within tight specifications, especially for an optical
telescope like the TNG. However, a good design of the slewing phase is
essential as well. The slew phase must be as short as possible, so ideally the
telescope axes should be driven at maximum speed, but this is impossible
because of physical limits for the next position derivatives, i.e.
acceleration and jerk. Overall, a system approach for servo control shall be
adopted, addressing simultaneously all issues related to both slew and
tracking phases.

With reference to the PI (Proportional-Integral) controller, one of the most
popular tracking position controllers for telescopes as well, one problem is
represented by the large overshoot caused by the integral action during the
slew phase, when the position error is large. One workaround is the
implementation of different position controllers for slew and tracking, with a
transition phase between the two; this was the first solution implemented at
the VLT as described in Ref. 9 (later the implementation was changed as
described in Ref. 10 adopting a command preprocessor, see also Section 5.1).

A more sophisticated and flexible solution, still attractively easy, is the
implementation of a variable structure controller in the loop (fig. 6),
pioneered in the telescope field at the TNG itself in the nineties. With a
variable structure controller, the parameters are dynamically changed as a
function of the position error, with the following guidelines.

  * •

When the axis is tracking and the error is small, the controller must
guarantee optimal performance and minimal servo error.

  * •

When the error is large, using the tracking gains would cause instability. The
position gain must be smaller and the integral gain must be zero for errors
larger than a threshold.

  * •

Between the two extreme conditions (large and small errors), it shall be
possible to shape the gain functions as appropriate to allow for a smooth
transition.

The implementation of these principles enables a smooth switch from a purely
proportional to a proportional-integral controller, varying continuously the
proportional and integral gains over the transition phase.

In the TNG case, this has been implemented through a simple linear
interpolation of three points in the gain-error plane; the principle is simply
described by fig. 7. The variable gain function is fully described by the
coordinates of the three points, i.e. by 6 parameters. Thus, the PI variable
structure controller is described by 12 parameters, fully defining Kp(e) and
Ki(e). With reference to the transition phase, it is worth mentioning that
more mathematically elegant connections of the two extreme points would not
give, for the shaping of the functions, the same total freedom of a simple
linear interpolation of points. In the TNG case, setting only one intermediate
point has been enough.

---
Figure 6:  Position loop adopting a Variable Structure controller.
---
Figure 7:  Variable structure controller gain. The K(e) function comes from
the linear interpolation of three points. The two extremes define the gain for
the slewing phase (large errors) and the tracking phase (small errors); the
central point gives the necessary degrees of freedom for the transition phase.
---
Figure 8:  Position loop adopting a command preprocessor in input.

##  5 Alternative implementations

###  5.1 Command preprocessor

Other approaches to the problems outlined in section 4 are possible, rather
than recurring to variable structure controllers. An effective solution can be
designing the slew phase by pre-filtering the axes reference positions coming
from the astrometric sequence, a solution sometimes referred to as command
pre-preprocessor or command shaper (Fig. 8). This way, the designer can
constraint the system variables by planning a trajectory that will never
exceed the velocity and acceleration limits, with two consequences: the
telescope dynamics will always work in the linear regime, as any saturation
will be prevented by design; the position error will always be limited,
effectively preventing the overshoot problem caused by the integral action.

The trajectory planning is not trivial because the telescope target is moving,
although slowly, so the final destination of the telescope slew is not
rigorously known in advance. Of course, in theory the whole trajectory of the
celestial object might be precomputed; but this would require the servo
control software embeds or replicates the astrometry computations.

This is normally an impractical solution, because the two tasks of astrometry
and servo control are usually decoupled and implemented separately: this is
the case at the TNG too. An astrometry module computes the axes coordinates
and generates the set-points for the servo control; on the servo side, the
astrometry is seen as a black-box that outputs periodically the coordinates to
track. Thus, the instantaneous trajectory set-points are known, but all the
future ones are unknowns. Therefore, an essential feature for a command
preprocessor at the TNG is the ability to work recursively, having no
knowledge of the future positions but using instantaneous inputs only.

Examples of slew trajectory planners are in Ref. 11, 12, 13, 14, 15, 16, 17.
They differ significantly in the approach and implementation. Here, the
preference for a possible future implementation at the TNG goes for a
recursive scheme guaranteeing the continuity of velocity and acceleration,
i.e. for a preprocessor implementing a C2 class trajectory.

The implementation of a command shaper may lead to the implementation of a
single position controller for any phase (large, medium and small errors),
removing the need of variable gains.

###  5.2 Model-based control

The current PI controller improves the performance when the gains are
increased, but makes the system unstable if the gains are too high, causing
vibrations. The reason is that the telescope is not totally rigid but has
flexible modes. A controller that suppresses the telescope vibrations when the
gains are increased would be helpful. This is what a model-based controller
aims to do. Therefore, although the current implementation fulfills the
requirements, moving to a model-based controller could in principle further
enhance the system.

One option is the LQG control. In the optimal control theory, a LQG controller
is a combination of a Kalman filter (a Linear Quadratic state Estimator - LQE)
with a Linear Quadratic Regulator (LQR). The tracking version of the LQG
controller described in Ref. 7 has some similarities with the PI control
scheme and may be a candidate for an experimentation at the TNG.

As shown in Fig. 9, the LQG is implemented in the position loop and the
”Telescope Axis”, to be modelled through a white-noise identification
procedure, is the velocity closed-loop. Here, the command preprocessor is
mandatory: in fact the capability to shape the trajectory not exceeding
velocity and acceleration limits allows to work always in the linear regime,
as needed by this method. The scheme includes an estimator of the states of
the flexible modes from the axis position and the speed command (output of
position controller). The flexible mode controller $K_{f}$ action is added to
the PI output, allowing the suppression of the telescope vibrations. The
tracking version of the LQG controller is a PI augmented with a flexible mode
controller, where the increase of the proportional gain cause no longer
instability because the vibrations are suppressed.

---
Figure 9:  Position loop adopting a command preprocessor and a LQG controller.

##  6 Conclusions

Within the framework of a hardware modernization effort at the TNG, more than
two decades into the start of telescope operations, the servo control
algorithms are consequently under update, to cope with the new drive system.
After the control algorithms redesign and the validation of the tracking
performance, the azimuth axis completed the transition to the new hardware.
The open-loop axis model was identified and the servo controllers were
redesigned from scratch. For the servo control algorithm redesign, few days of
technical time in day-hours were allocated at the observatory, without
stopping the telescope night operations.

The minimum goal was leaving the telescope operational with the new hardware
system, with at least the same satisfactory performance shown for the previous
twenty years. This goal was achieved, as the final performance exceeded the
previous ones.

Given the limited time allocated for this first intervention, it was
conservatively decided to start from the control philosophy implemented by the
same authors at telescope commissioning time. New options are under study,
which might further improve the TNG tracking quality.

## References

  * [1] Gonzalez, M., Ventura, H. P., Diaz, M. H., Schipani, P., Ghedina, A., Gómez, J. S. J., Hernandez, N., Guerra, J., Lodi, M., Rodríguez, A. L. R., Riverol, C., Gonzalez, C., Galli, A., Colapietro, M., and Savarese, S., “Telescopio Nazionale Galileo control system upgrade, new milestone achieved: Azimuth axis completed,” Proc. SPIE 11445 (2020).
  * [2] Abbott, T. M. C., Schumacher, G., Warner, M., Mondaca, E., Schmidt, R. E., Cantarutti, R., Martinez, M., Estay, O., Delgado, F., and Walker, A., “Cerro Tololo Inter-American observatory, Victor M. Blanco 4-m telescope: an upgrade to the telescope control system,” Proc. SPIE 7012, 70123K (2008).
  * [3] Sprayberry, D., Dunlop, P., M.Evatt, Reddell, L., Gott, S., George, J. R., Donaldson, J., Stupak, R. J., R.Marshall, Abareshi, B., Stover, D., Warner, M., Cantarutti, R. E., and Probst, R. G., “Modernization of the Mayall telescope control system: design, implementation, and performance,” Proc. SPIE 9906, 990608 (2016).
  * [4] Krasuski, T., Thompson, P., McCarney, B., Kwok, S., James, E., Tsubota, K., Pollard, M., and McCann, K., “Implementation of W. M. Keck observatory’s telescope control system upgrade,” Proc. SPIE 10700, 107005I (2018).
  * [5] Erm, T. and Gutierrez, P., “Integration and tuning of the VLT drive systems,” Proc. SPIE 4004, 490–499 (2000).
  * [6] Schipani, P., Arcidiacono, C., , Argomedo, J., Dall’Ora, M., D’Orsi, S., Farinato, J., Magrin, D., Marty, L., Ragazzoni, R., and Umbriaco, G., “The tracking control system of the VLT Survey Telescope,” Rev. Sci. Instrum. 83, 094501 (2012).
  * [7] Gawronski, W., [Modeling and control of antennas and telescopes ], Springer (2008).
  * [8] Gawronski, W., [Adavanced structural dynamics and active control of structures ], Springer (2004).
  * [9] Chiesa, M. and Sandrock, S., “VLT-TRE-ESO-17130-1211,” (2001).
  * [10] Sandrock, S., Di Lieto, N., Pettazzi, L., and Erm, T., “Design and implementation of a general main axis controller for the eso telescopes,” Proc. SPIE 8451, 84510N (2012).
  * [11] Tyler, S. R., “A trajectory preprocessor for antenna pointing,” TDA Progress Report 42-118, 139–159 (1994).
  * [12] Olberg, M., Lindeborg, C., Seyf, A., and Kastengren, C. F., “Simple, robust digital controller for the onsala 20-m radio telescope,” Proc. SPIE 2479, 257–265 (1995).
  * [13] Prestage, R. M., Oliveira, F. J., Bailey, J. A., Friberg, P., and Kenderdine, S., “Route-finding algorithm for the optimal control of telescope slews,” Proc. SPIE 3112, 76–87 (1997).
  * [14] Gawronski, W. and Almassy, W. T., “Command preprocessor for radio telescopes and microwave antennas,” IEEE Antennas Propag. Mag. 44 (2), 30–37 (2002).
  * [15] Smith, D. R., “A polynomial-based trajectory generator for improved telescope control,” Proc. SPIE 7019, 701909 (2008).
  * [16] Hunter, T. R., Wilson, R. W., Kimberk, R., Leiker, P. S., Patel, N. A., Blundell, R., Christensen, R. D., Diven, A. R., Maute, J., Plante, R. J., Riddle, P., and Young, K. H., “The digital motion control system for the submillimeter array antennas,” J. Astron. Instrum. 2 (1), 1350002 (2013).
  * [17] Savarese, S., Perrotta, F., Schipani, P., Schillirò, F., Capasso, G., Colapietro, M., and D’Orsi, S., “Trajectory generation methods for radio and optical telescopes,” Proc. SPIE 11445 (2020).


# Jacobi multipliers and Hojman symmetry

José F. Cariñena and Manuel F. Rañada
Departamento de Física Teórica and IUMA, Universidad de Zaragoza.
50009, Zaragoza, Spain
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

The geometric intrinsic approach to Hojman symmetry is developed and use is
made of the theory of the Jacobi last multipliers to find the corresponding
conserved quantity for non divergence-free vector fields. The particular cases
of autonomous Lagrangian and Hamiltonian systems are studied as well as the
generalization of these results to normalizer vector fields of the dynamics.
The nonautonomous cases, where normalizer vector fields play a relevant role,
are also developed.

MSC: 70G45, 70H03,70H05.

PACS: 02.30.Hq, 02.40.-k, 45.20.Jj, 76M60, 87.10.Ed.

Keywords: Hojman symmetry, Jacobi multipliers

##  1 Introduction

The first-integrals of a vector field in a $n$-dimensional manifold $M$,
$X={\displaystyle\sum_{i=1}^{n}X^{i}(x)\partial/\partial x^{i}}$, play a very
relevant role in its integrability, and therefore in that of the corresponding
autonomous system of first-order differential equations:

$\dot{x}^{i}=X^{i}(x),\quad i=1,\ldots,n.$ (1.1)

Moreover, the knowledge of other tensor fields that are invariant under the
given vector field furnishes tools for reduction to simpler cases, and the
best known example is the case of an invariant symplectic form $\omega$, for
which one of the more celebrated related results is that of the Noether
theorem, that in this case amounts to say that if a Hamiltonian vector field
$X_{F}$ is a symmetry of the Hamiltonian function $H$, then the function $F$
is a constant of the motion for the dynamical Hamiltonian vector field
$X_{H}$. The very simple proof is based on the skew-symmetry of the Poisson
bracket defined by the symplectic structure, because
$X_{H}F=\\{F,H\\}=-X_{F}H$. The particular case of regular Lagrangians can be
considered as a particular case of a Hamiltonian system where, in this case,
the symplectic form $\omega_{L}$ and the Hamiltonian function $H$, the energy
function, are defined from the Lagrangian $L$. The symmetries of both the
symplectic form and the energy function correspond to the invariance of the
Lagrangian, may be up to addition of a gauge term. Just to mention another
instance, the existence of an invariant volume form for a vector field $X$ in
an oriented manifold, which is equivalent to know a Jacobi multiplier for $X$,
was used by Jacobi to prove that if $(n-2)$ first integrals are known, the
vector field is integrable by quadratures.

Even if the usual way for finding first-integrals is the method suggested by
Noether theorem of introducing symmetries, there is an alternative method for
finding non-Noether constants of motion that is not based on symmetries but on
the existence of alternative compatible geometric structures for the
description of the vector field, what leads to the existence of a recursion
operator. The case of the Lagrangian formalism for a system with a one-
dimensional configuration space was given in [1], while the multidimensional
case was given in [2] and the geometric approach in intrinsic terms was given
in [3].

The third approach we want to study here is the one started in a coordinate
dependent way by Hojman [4] for the case of a divergence-free non-autonomous
system of second-order ordinary differential equations and later on extended
by González-Gascón in [5] for a non divergence-free case, and that it is
becoming more and more important during the last years for its applications in
$f(R)$-gravity and FRW cosmology [6, 7, 8, 9, 10, 11, 12]. The geometric
approach to autonomous systems provides a method to incorporate holonomic
constraints by replacing the systems by vector fields in a differentiable
manifold and allows to develop coordinate free formulations of the problem. In
this way the geometric theory of reduction and symmetry may be very helpful.
Therefore, our aim in this paper is to present a geometric approach
generalising the previously known cases and unifying them in a common simple
geometric formulation, by showing that they appear as particular cases of a
more general theory.

The paper is organized as follows: Section 2 summarises basic results of the
differential geometry concerning the concept of divergence of a vector field
and a basic relation between the divergences of two vector fields, and we also
recall the fundamental concept of Jacobi multiplier introduced by him to
integrate by quadratures a vector field when $(n-2)$ functionally independent
first-integrals and a Jacobi multiplier are known. Using these concepts we can
introduce some mathematical results which can be used to obtain a first-
integral of the vector field $X$ in case that an infinitesimal symmetry $Y$ of
$X$ or a normalizer vector field are known. After this geometric presentation
we particularise in Section 3 to the cases of autonomous Lagrangian and
Hamiltonian systems, while nonautonomous systems of first-order differential
equations are analysed in Section 4, and in particular the Hamilton systems.
Section 5 is devoted to nonautonomous systems of second-order differential
equations with a special emphasis on those admitting a regular Lagrangian
formulation. Many results on constants of motion scattered by the physics
literature arise here as particular cases of this more general unifying
presentation.

##  2 First integrals derived from infinitesimal symmetries of vector fields

Different relations among vector fields and their divergences can be used to
establish first integrals and integral invariants for vector fields in a
manifold $M$. Given a vector field $X\in\mathfrak{X}(M)$ in an oriented
manifold $(M,\Omega)$, we denote ${\rm div}(X)$ the function defined by
$\mathcal{L}_{X}\Omega={\rm div}(X)\,\Omega$. In the particular case of a
choice of local coordinates $(x^{1},\ldots,x^{n})$ such that
$\Omega=dx^{1}\wedge\cdots\wedge dx^{n}$, if the local expression of $X$ is
$X={\displaystyle\sum_{i=1}^{n}X^{i}\frac{\partial}{\partial x^{i}}}$, then
the expression of ${\rm div}(X)$ is the usual one:

${\rm div}(X)=\sum_{i=1}^{n}\frac{\partial X^{i}}{\partial x^{i}}.$

Note that for each function $f\in C^{\infty}(M)$, from
$\mathcal{L}_{f\,X}\Omega=\mathcal{L}_{X}(f\,\Omega)=(X(f)+f\,{\rm
div}(X))\Omega$ we have that ${\rm div}(fX)=X(f)+f\,{\rm div}(X)$.

There is a geometric relation playing a fundamental rôle: If $X,Y$ is an
arbitrary pair of vector fields in an oriented manifold $(M,\Omega)$, then

$\mathcal{L}_{X}({\rm div\,}(Y))-\mathcal{L}_{Y}({\rm div\,}(X))={\rm
div\,}([X,Y]),$ (2.1)

because as the Lie derivatives of a $p$-form $\alpha$ satisfy the relation

$\mathcal{L}_{X}(\mathcal{L}_{Y}\alpha)-\mathcal{L}_{Y}(\mathcal{L}_{X}\alpha)=\mathcal{L}_{[X,Y]}\alpha,$
(2.2)

in the particular case $\alpha=\Omega$, we have

$\mathcal{L}_{X}(\mathcal{L}_{Y}\Omega)-\mathcal{L}_{Y}(\mathcal{L}_{X}\Omega)=\mathcal{L}_{[X,Y]}\Omega,$
(2.3)

and from the difference of

$\mathcal{L}_{X}(\mathcal{L}_{Y}\Omega)=\mathcal{L}_{X}({\rm
div}(Y)\,\Omega)=\mathcal{L}_{X}({\rm div}(Y))\,\Omega+{\rm div}(Y)\,{\rm
div}(X)\,\Omega$

and the corresponding relation

$\mathcal{L}_{Y}(\mathcal{L}_{X}\Omega)=\mathcal{L}_{Y}({\rm
div}(X)\,\Omega)=\mathcal{L}_{Y}({\rm div}(X))\,\Omega+{\rm div}(X)\,{\rm
div}(Y)\,\Omega$

we obtain (2.1).

The following theorem makes use of this geometric property and generalises
some previous results which had been obtained in [4, 13] by making use of a
coordinate dependent approach:

###### Theorem 1

Let $X$ be a divergence-free vector field and $Y$ an infinitesimal symmetry of
$X$. Then ${\rm div\,}(Y)$ is a constant of the motion for $X$.

Proof.- Let us apply the relation (2.1) to a divergence-free vector field $X$.
Then if the vector field $Y\in\mathfrak{X}(M)$ is a symmetry of $X$, it
satisfies $[Y,X]=0$ and therefore we obtain that $\mathcal{L}_{X}({\rm
div}(Y))=0$.

$\Box$

Next, in the following two theorems, we present two generalizations of this
result. First, we delete the condition for $X$ to be a divergence-free vector
field and second, the condition of $Y$ being an exact symmetry is substituted
for the weaker condition of symmetry of the one-dimensional distribution
generated by $X$. So, if the condition for $X$ to be a divergence-free vector
field is eliminated we can use the Jacobi multiplier theory and then we have a
more general result obtained in [5, 14]. Recall first that a positive function
$R\in C^{\infty}(M)$ is a Jacobi multiplier for $X\in\mathfrak{X}(M)$ in the
oriented manifold $(M,\Omega)$ if the vector field $R\,X$ is such that
$\mathcal{L}_{R\,X}\Omega=0$, or equivalently, as
$\mathcal{L}_{R\,X}(\Omega)=d(i(R\,X)\Omega=d(i(X)(R\,\Omega))=\mathcal{L}_{X}(R\,\Omega)$,
if $\mathcal{L}_{X}(R\,\Omega)=0$. Note that

$\mathcal{L}_{X}(R\,\Omega)=\mathcal{L}_{X}(R)\,\Omega+R\,\mathcal{L}_{X}\Omega=(\mathcal{L}_{X}(R)+R\,{\rm
div}(X))\Omega,$

and therefore the condition for $R$ to be a Jacobi multiplier (see [15] and
references therein for more details) becomes

$\mathcal{L}_{X}(R)+R\,{\rm div}(X)=0\Longleftrightarrow{\rm
div}(X)+\mathcal{L}_{X}(\log R)=0.$ (2.4)

Remark that from the above mentioned property
$\mathcal{L}_{X}(R\,\Omega)=\mathcal{L}_{R\,X}\Omega$, we see that $R$ is a
Jacobi multiplier for $X$ if and only if the vector field $X$ is divergence-
free with respect to $R\,\Omega$, i.e., we see the equivalence of searching
for Jacobi multipliers on one side and the more geometric concept of invariant
volume-forms on the other: a volume form $\Omega^{\prime}=R\,\Omega$ is
$X$-invariant, i.e. $\mathcal{L}_{X}\Omega^{\prime}=0$, if and only if $R$ is
a Jacobi multiplier for the vector field $X$ in the oriented manifold
$(M,\Omega)$. Moreover, this property also shows that $R$ is a Jacobi
multiplier for $X$ in the oriented manifold $(M,\Omega)$ iff and only if
$f\,R$ is a Jacobi multiplier for $X$ in the oriented manifold
$(M,f^{-1}\,\Omega)$, for each positive function $f$.

###### Theorem 2

If the function $R$ is a Jacobi multiplier for the vector field $X$ and the
vector field $Y$ is an infinitesimal symmetry of $X$, i.e. $[X,Y]=0$, then the
function

$I={\rm div\,}(Y)+Y(\log R)$ (2.5)

is a constant of the motion for $X$.

Proof.- If $[X,Y]=0$, then according to the relation (2.1) for such a pair
$X,Y,$ of vector fields we have that $\mathcal{L}_{X}{\rm
div\,}(Y)=\mathcal{L}_{Y}{\rm div\,}(X)$. If $R$ is a Jacobi multiplier for
the vector field $X$, then (2.4) shows that
$\textrm{div}(X)=-\mathcal{L}_{X}(\log R)$, and consequently in this case the
equation (2.1) can be rewritten as

$\mathcal{L}_{X}{\rm div\,}(Y)=-\mathcal{L}_{Y}(\mathcal{L}_{X}(\log
R))=-\mathcal{L}_{X}(\mathcal{L}_{Y}(\log R)).$

Therefore, $\mathcal{L}_{X}I=0$.

$\Box$

Another important case is when the vector field $Y$ is not a symmetry of the
vector field $X$ but only of the one-dimensional distribution spanned by $X$,
i.e. there exists a function $h\in C^{\infty}(M)$ such that $[Y,X]=h\,X$. In
this case the results of the previous theorems can be generalised to:

###### Theorem 3

i) If the vector field $X$ is divergence-free and the vector field $Y$ is an
infinitesimal symmetry of the one-dimensional distribution generated by $X$,
i.e. $[Y,X]=h\,X$, then the function $\textrm{div\,}(Y)+h$ is a constant of
the motion for $X$.

ii) If the function $R$ is a Jacobi multiplier for the vector field $X$ and
the vector field $Y$ is an infinitesimal symmetry of the one-dimensional
distribution generated by $X$, i.e. i.e. there exists a function $h\in
C^{\infty}(M)$ such that $[Y,X]=h\,X$, then the function

$I={\rm div\,}(Y)+Y(\log R)+h$ (2.6)

is a constant of the motion for $X$.

Proof.- i) If ${\rm div\,}(X)=0$, then $\mathcal{L}_{X}\Omega=0$, and hence,
using (2.1),

$\mathcal{L}_{X}(\mathcal{L}_{Y}\Omega)=\mathcal{L}_{Y}(\mathcal{L}_{X}\Omega)+\mathcal{L}_{-h\,X}\Omega=-\mathcal{L}_{X}(h\,\Omega)=-X(h)\,\Omega,$

and therefore, from

$\mathcal{L}_{X}(\mathcal{L}_{Y}\Omega)=\mathcal{L}_{X}({\rm
div\,}(Y)\,\Omega)=\mathcal{L}_{X}({\rm div\,}(Y))\Omega$

we obtain:

$\mathcal{L}_{X}({\rm div\,}(Y)+h)=0.$

ii) Using the relation (2.1) for the vector fields $X$ and $Y$, and the Jacobi
multiplier condition ${\rm div\,}(X)+X(\log R)=0$, together with $[Y,X]=h\,X$,
we find that

$\mathcal{L}_{X}{\rm div\,}(Y)-\mathcal{L}_{Y}(-\mathcal{L}_{X}(\log R))=-{\rm
div\,}(h\,X)=-X(h)-h\,{\rm div\,}(X),$

and as the commutator $[\mathcal{L}_{X},\mathcal{L}_{Y}]$ is
$\mathcal{L}_{-h\,Y}$, the preceding relation can be rewritten as

$\mathcal{L}_{X}({\rm div\,}(Y)+Y(\log R))+\mathcal{L}_{h\,X}(\log
R)=-X(h)+h\,X(\log R)$

and simplifying terms,

$\mathcal{L}_{X}({\rm div\,}(Y)+Y(\log R)+h)=0.$

$\Box$

Given a vector field $X\in\mathfrak{X}(M)$ the vector fields
$Y\in\mathfrak{X}(M)$ that are symmetries of the distribution generated by $X$
are sometimes called normalizers of $X$.

With the particular choice of local coordinates $(x^{1},\ldots,x^{n})$ such
that $\Omega=dx^{1}\wedge\cdots\wedge dx^{n}$, if the local expressions of the
vector fields $X$ and $Y$ are

$X=\sum_{i=1}^{n}X^{i}\frac{\partial}{\partial x^{i}},\qquad
Y=\sum_{i=1}^{n}Y^{i}\frac{\partial}{\partial x^{i}},$

the coordinate expression of constant of motion (2.6) is

$I=\sum_{i=1}^{n}\frac{\partial Y^{i}}{\partial
x^{i}}+\sum_{i=1}^{n}Y^{i}\frac{\partial(\log R)}{\partial x^{i}}+h\,,$

and the condition $[Y,X]=h\,X$ means that

$Y(X^{i})-X(Y^{i})=h\,X^{i}.$

##  3 Hojman symmetry in the Lagrangian and Hamiltonian formalism

We can extend the preceding results to systems of second-order differential
equations in normal form,

$\ddot{x}^{i}=F^{i}(x,\dot{x}),\quad i=1,\ldots,n,$ (3.1)

as done in [5, 14], and the particular case in which this system is equivalent
to the set of Euler-Lagrange equations of a regular Lagrangian is of a special
interest.

Remark that such a system of $n$ second-order differential equation can be
associated with a system of $2n$ first-order differential equations

$\left\\{\begin{array}[]{rcl}\dot{x}^{i}&=&v^{i}\\\
\dot{v}^{i}&=&F^{i}(x,v)\end{array}\right.,\quad i=1,\ldots,n,$ (3.2)

that as a particular case of (1.1) has associated the vector field

$\Gamma=\sum_{i=1}^{n}\left(v^{i}\frac{\partial}{\partial
x^{i}}+F^{i}(x,v)\frac{\partial}{\partial v^{i}}\right).$ (3.3)

When the configuration space $Q$ of the system is a $n$-dimensional Euclidean
space the Euclidean coordinates $(x^{1},\ldots,x^{n})$ induce global
coordinates in its tangent bundle, denoted
$(x^{1},\ldots,x^{n},v^{1},\ldots,v^{n})$ and there is a natural volume
element (tangent and cotangent bundles can be identified by the Euclidean
metric) given by $\Omega=dx^{1}\wedge\cdots\wedge dx^{n}\wedge
dv^{1}\wedge\cdots\wedge dv^{n}$. If the dynamical vector field (3.3) is
determined by the Lagrangian $L$, i.e. $i(\Gamma)\omega_{L}=dE_{L}$ (see e.g
[16]), we know (see e.g. [15] and references therein) that a particular Jacobi
multiplier for given by the determinant of the Hessian matrix in the
velocities, with elements $W_{ij}=\partial^{2}L/{\partial v^{i}\partial
v^{j}}$. Therefore, the constant of motion obtained in [14] is just a
particular case of the expression (2.5) for $R$ equal to the determinant of
the Hessian matrix $W$. In the more general case of $Q$ a $n$-dimensional
manifold, a local chart of coordinates $(x^{1},\ldots,x^{n})$ for $Q$ induces
a local chart of coordinates in its tangent bundle, denoted
$(x^{1},\ldots,x^{n},v^{1},\ldots,v^{n})$, and an associated volume element in
such a chart given by $\Omega=dx^{1}\wedge\cdots\wedge dx^{n}\wedge
dv^{1}\wedge\cdots\wedge dv^{n}$. If the dynamical vector field $\Gamma$ given
in such coordinates by (3.3) is determined by the Lagrangian $L$, i.e.
$i(\Gamma)\omega_{L}=dE_{L}$, which implies $\mathcal{L}_{\Gamma}\omega_{L}=0$
(see e.g [16]), we know (see e.g. [15] and references therein) that a
particular Jacobi multiplier for $\Gamma$ with respect to the volume form
$\Omega$ is given by the determinant of the Hessian matrix $W$ in the
velocities, with elements $W_{ij}=\partial^{2}L/{\partial v^{i}\partial
v^{j}}$. Actually, from $\mathcal{L}_{\Gamma}\omega_{L}=0$, we see that
$(\omega_{L})^{\wedge n}$ is an invariant volume under $\Gamma$, and the
proportionality factor of $(\omega_{L})^{\wedge n}$ and $\Omega$ is $\det W$.
Therefore, the constant of motion obtained in [14] is just a particular case
of the expression (2.5) for $R$ equal to the determinant of the Hessian matrix
$W$.

Let us now consider the coordinate expression of the vector field $Y$ defined
in $TQ$ with local coordinates $(x^{i},v^{i})$

$Y=\sum_{i=1}^{n}\left(Y^{i}(x,v)\frac{\partial}{\partial
x^{i}}+{\bar{Y}^{i}}(x,v)\frac{\partial}{\partial v^{i}}\right),$ (3.4)

where the fuctions $\bar{Y}^{i}$ are given by $\bar{Y}^{i}=\Gamma(Y^{i})$,
because of the symmetry condition $[\Gamma,Y]=0$. Then the coordinate
expressions of the constants of motion $I$ given in Theorem 1 and in (2.5) of
Theorem 2 for the vector field $\Gamma$ are given, respectively, by

$I=\sum_{i=1}^{n}\left(\frac{\partial Y^{i}}{\partial
x^{i}}+\frac{\partial\bar{Y}^{i}}{\partial v^{i}}\right),$ (3.5)

and

$I=\sum_{i=1}^{n}\left(\frac{\partial Y^{i}}{\partial
x^{i}}+\frac{\partial\bar{Y}^{i}}{\partial
v^{i}}\right)+\frac{1}{R}\sum_{i=1}^{n}\left(Y^{i}\frac{\partial R}{\partial
x^{i}}+{\bar{Y}^{i}}\frac{\partial R}{\partial v^{i}}\right),$ (3.6)

These expressions coincide with those of [6, 7].

Similarly, let us now consider the time evolution in phase space given in
local coordinates $(q^{1},\cdots,q^{n},p_{1},\cdots,p_{n})$ of the cotangent
bundle $\pi:T^{*}Q\to Q$ induced from a chart $(q^{1},\ldots,q^{n})$ in the
base $Q$ by

$\left\\{\begin{array}[]{rcl}\dot{q}^{i}&=&Q^{i}(q,p)\\\
\dot{p}_{i}&=&P_{i}(q,p)\end{array}\right.,$ (3.7)

whose solutions correspond to the integral curves of the vector field

$X=\sum_{i=1}^{n}\left(Q^{i}(q,p)\,\frac{\partial}{\partial
q^{i}}+P_{i}(q,p)\,\frac{\partial}{\partial p_{i}}\right)\,.$

and let $\Omega$ be the volume form induced from exterior product of its
natural symplectic structure, $\Omega=dq^{1}\wedge\cdots dq^{n}\wedge
dp_{1}\wedge\dots dp_{n}$. Recall that the cotangent bundle is the local
prototype of an exact symplectic manifold. There is a canonical 1-form
$\theta\in\bigwedge^{1}(T^{*}Q)$ defined by
$\theta_{\alpha}=\pi^{*}_{\alpha}\alpha$, with local expression
$\theta={\displaystyle\sum_{i=1}^{n}}p_{i}\,dq^{i}$. Then $\omega=-d\theta$ is
a symplectic form on $T^{*}Q$ and the local coordinates are Darboux
coordinates.

Then if $R$ is a Jacobi multiplier for $X$, satisfying therefore the equation
(2.4),

$\sum_{i=1}^{n}\left(\frac{\partial Q^{i}}{\partial q^{i}}+\frac{\partial
P_{i}}{\partial p_{i}}\right)+X(\log R)=0,$

and the vector field $Y$, given in these local coordinates by

$Y={\displaystyle\sum_{i=1}^{n}}\left(\xi^{i}(q,p)\frac{\partial}{\partial
q^{i}}+\eta_{i}(q,p)\,\frac{\partial}{\partial p_{i}}\right),$

is a symmetry of the vector field $X$, the constant of motion (2.5) for $X$ is
in this case

$I=\sum_{i=1}^{n}\left(\frac{\partial\xi^{i}}{\partial
q^{i}}+\frac{\partial\eta_{i}}{\partial p_{i}}+\xi^{i}\frac{\partial(\log
R)}{\partial q^{i}}+\eta_{i}\frac{\partial(\log R)}{\partial p_{i}}\right),$

as found in [17]. Of course, if $(M,\omega)$ is a symplectic manifold and $X$
is a Hamiltonian vector field and we consider the volume form in $M$ obtained
by exterior products of the symplectic form, as $X$ is divergence-free, the
multiplier $R$ is just a constant.

Finally, if the vector field $Y$ is a normalizer for $X$, i.e., there exists a
fuction $h\in C^{\infty}(M)$ such that $[Y,X]=h\,X$, then the constant of
motion corresponding to that of Theorem 3 is

$I=\sum_{i=1}^{n}\left(\frac{\partial\xi^{i}}{\partial
q^{i}}+\frac{\partial\eta_{i}}{\partial p_{i}}+\xi^{i}\frac{\partial(\log
R)}{\partial q^{i}}+\eta_{i}\frac{\partial(\log R)}{\partial p_{i}}\right)+h.$

##  4 Hojman symmetry for nonautononous systems of first-order differential
equations

In the geometric approach to non-autonomous systems of first-order
differential equations,

$\dot{x}^{i}=X^{i}(t,x),\quad i=1,\ldots,n,$ (4.1)

time is a new coordinate and then we can change the parameter of the solution
curves. This amounts to consider the curves not as parametrized curves but
simply as one-dimensional submanifolds and to consider as equivalent two
curves that are related by a reparametrization. Therefore we associate with
such a non-autonomous system a one-dimensional, and therefore integrable,
distribution, generated either by a vector field $X$ in $\mathbb{R}\times M$,
or by any other proportional one, $f\,X$, with $f$ a nonvanishing function, in
such a way that its one-dimensional integral manifolds are the equivalence
class of integral curves of a representative of the vector field. We can
choose, for instance

$X=\frac{\partial}{\partial
t}+\sum_{i=1}^{n}X^{i}(t,x)\frac{\partial}{\partial x^{i}},$ (4.2)

if for a vector field of the distribution we have $\langle dt,X\rangle\neq 0$.

The concept of infinitesimal symmetry must be changed and we say that the
vector field $Y\in\mathfrak{X}(\mathbb{R}\times M)$ is a symmetry when the
image under its local flow of an integral curve of $X$ is, up to a
reparametrization, an integral curve of $X$. This condition means that there
exists a function $\lambda\in C^{\infty}(\mathbb{R}\times M)$ such that
$[Y,X]=\lambda\,X$, or in other terms, that the one-dimensional distribution
generated by $X$ is invariant under the flow of the vector field $Y$. Note
that if the coordinate expression of the vector field $Y$ is

$Y=Y^{0}(t,x)\frac{\partial}{\partial
t}+\sum_{i=1}^{n}Y^{i}(t,x)\frac{\partial}{\partial x^{i}},$

then the condition $[Y,X]=\lambda\,X$ implies that $\lambda=-X(Y^{0})$, as a
consequence of the explicit form (4.2) of the representative vector field $X$.
In this non-autonomous case, if either $M$ is an Euclidean space or choosing a
local chart in $M$, we have the associated the volume form in
$\mathbb{R}\times M$ given by $\Omega=dt\wedge dx^{1}\wedge\cdots\wedge
dx^{n}$ and then ${\rm div\,}(X)$ and ${\rm div\,}(Y)$ take the form in such
coordinates

${\rm div\,}(X)=\sum_{i=1}^{n}\frac{\partial X^{i}}{\partial
x^{i}},\qquad\qquad{\rm div\,}(Y)=\frac{\partial Y^{0}}{\partial
t}+\sum_{i=1}^{n}\frac{\partial Y^{i}}{\partial x^{i}}.$

The results corresponding to Theorem 1, having in mind the new concept of
symmetry is, according to i) of Theorem 3:

###### Theorem 4

i) If the vector field $X$ is divergence-free and the vector field $Y$ is an
infinitesimal symmetry of the one-dimensional distribution generated by $X$,
then the function $\textrm{div\,}(Y)-X(Y^{0})$ is a constant of the motion for
$X$.

ii) If the function $R$ is a Jacobi multiplier for the vector field $X$ and
the vector field $Y$ is an infinitesimal symmetry of the one-dimensional
distribution generated by $X$, then the function

$I={\rm div\,}(Y)+Y(\log R)-X(Y^{0})$ (4.3)

is a constant of the motion for $X$.

$\Box$

In the particular case of $M=T^{*}Q$, if we choose a chart for
$\mathbb{R}\times T^{*}Q$ induced from a local chart for $Q$ and the local
expression of the vector field $X$ is

$X=\frac{\partial}{\partial
t}+\sum_{i=1}^{n}\left(Q^{i}(t,q,p)\frac{\partial}{\partial
q^{i}}+P_{i}(t,q,p)\frac{\partial}{\partial p_{i}}\right),$

which corresponds to a non-autonomous system as (3.7) but where now the
functions $Q$ and $P$ also depend on $t$, if the coordinate expression of a
vector field $Y$ on $\mathbb{R}\times T^{*}Q$ is

$Y=\sigma\,\frac{\partial}{\partial
t}+\sum_{i=1}^{n}\left(\xi^{i}\,\frac{\partial}{\partial
q^{i}}+\eta_{i}\,\frac{\partial}{\partial p_{i}}\right),$

then, as

$[Y,X]=-(X(\sigma)\,\frac{\partial}{\partial
t}+\sum_{i=1}^{n}\left((Y(Q^{i})-X(\xi^{i}))\frac{\partial}{\partial
q^{i}}+(Y(P_{i})-X(\eta_{i}))\frac{\partial}{\partial p_{i}}\right),$

the conditions for $Y$ to be a normalizer of $X$, $[Y,X]=h\,X$, are
$h=-X(\sigma)$ and $Y(Q^{i})-X(\xi^{i})=-X(\sigma)\,Q^{i}$ together with
$Y(P_{i})-X(\eta^{i})=-X(\sigma)\,P_{i}$, and in thus case the associated
constant of motion, for each Jacobi multiplier $R$ for $X$ with respect to the
volume form $dt\wedge dq^{1}\wedge\cdots\wedge dq^{n}\wedge dp_{1}\wedge\cdots
dp_{n}$ is

$I=\frac{1}{R}\frac{\partial(R\sigma)}{\partial
t}+\frac{1}{R}\sum_{i=1}^{n}\left(\frac{\partial(R\xi^{i})}{\partial
q^{i}}+\frac{\partial(R\eta_{i})}{\partial p_{i}}\right)-X(\sigma).$

what shows that the results indicated in [18] are but particular cases of the
more general theory developed here.

##  5 Hojman symmetry for nonautononous systems of second-order differential
equations

A non-autonomous systems of second-order differential equations,

$\frac{d^{2}x^{i}}{dt^{2}}=F^{i}(t,x,\dot{x}),\quad i=1,\ldots,n,$ (5.1)

has an associated system of $(2n+1)$ first-order differential equations

$\left\\{\begin{array}[]{rcl}{\displaystyle\frac{dt}{ds}}&=&1\\\
{\displaystyle\frac{dx^{i}}{ds}}&=&v^{i}\\\
{\displaystyle\frac{dv^{i}}{ds}}&=&F^{i}(t,x,v)\end{array}\right.\quad
i=1,\ldots,n,$ (5.2)

which is a particular case of (4.1) and it is therefore geometrically
described by a one-dimensional distribution $\mathcal{D}_{\Gamma}$ generated
by a vector field $\Gamma\in\mathfrak{X}(\mathbb{R}\times T\mathbb{R}^{n})$,
according to (4.2), given by

$\Gamma=\frac{\partial}{\partial
t}+\sum_{i=1}^{n}\left(v^{i}\frac{\partial}{\partial
x^{i}}+F^{i}(t,x,v)\frac{\partial}{\partial v^{i}}\right).$ (5.3)

Moreover, in the geometric approach $\ \mathbb{R}^{n}$ is replaced by a
$n$-dimensional manifold $Q$ and (5.3) is the local expression of a vector
field in the so-called evolution space $E=\ \mathbb{R}\times TQ$. Remark that
$\pi_{1}:\ \mathbb{R}\times Q\to\ \mathbb{R}$ defines a trivial bundle and the
jet-bundle of its sections is $J^{1}(\pi_{1})=\mathbb{R}\times TQ$.

This vector field $\Gamma$ can be characterized by the vanishing of the $2n$
1-forms in $\mathbb{R}\times TQ$

$\left\\{\begin{aligned} \theta^{i}&=dx^{i}-v^{i}\,dt,\\\
\Psi^{i}&=dv^{i}-F^{i}(t,x,v)\,dt,\end{aligned}\right.\quad i=1,\ldots,n,$
(5.4)

together with the condition $\langle dt,\Gamma\rangle=1$.

The first $n$ 1-forms $\theta^{i}$ are called contact 1-forms. Any other
representative of the distribution $\mathcal{D}_{\Gamma}$ will be of the form

$\bar{\Gamma}=\lambda(t,x,v)\Gamma=\lambda(t,x,v)\left(\frac{\partial}{\partial
t}+\sum_{i=1}^{n}\left(v^{i}\frac{\partial}{\partial
x^{i}}+F^{i}(t,x,v)\frac{\partial}{\partial v^{i}}\right)\right),$ (5.5)

with $\lambda(t,x,v)$ a nonvanishing function. Vector fields of type (5.5) are
called SODE vector fields while those of type (5.3) are said to be Newtonian
SODE (hereafter shortened as NSODE) vector fields.

We follow notation of [19] and we therefore refer to such a paper for
definitions. More details can be found in [20, 21]. An important point is that
the set of 1-forms

$\\{dt,\theta^{i}=dx^{i}-v^{i}\,dt,\Psi^{i}=dv^{i}-F^{i}(t,x,v)\,dt\\},$ (5.6)

is a local basis of 1-forms in $\mathbb{R}\times TQ$ which is the dual basis
of vector fields in $\mathfrak{X}(\mathbb{R}\times TQ)$, given by

$\left\\{\Gamma,\frac{\partial}{\partial x^{i}},\frac{\partial}{\partial
v^{i}}\right\\}.$ (5.7)

Note that while the contact forms only depend on the evolution space, the $n$
1-forms $\Psi^{i}$ depend also on the selected NSODE vector field $\Gamma$.
Particular examples for this case of second order differential equations were
studied in [4, 5, 13].

We only want to remark here that the Lagrangian formalism is defined for a
function $L\in C^{\infty}(\mathbb{R}\times TQ)$ by means of a $(1,1)$-tensor
in the evolution space $E$, called vertical endomorphism, whose local
coordinate expression is

$S=\sum_{i=1}^{n}\frac{\partial}{\partial
v^{i}}\otimes\theta^{i}=\sum_{i=1}^{n}\frac{\partial}{\partial
v^{i}}\otimes(dx^{i}-v^{i}\,dt).$ (5.8)

This may be used to define the 1-form $\Theta_{L}=L\,dt+dL\circ S$ in $E$,

$\Theta_{L}=L\,dt+\langle dL\circ
S,{\Gamma}\rangle\,dt+\sum_{i=1}^{n}\left\langle dL\circ
S,\frac{\partial}{\partial
x^{i}}\right\rangle\theta^{i}+\sum_{i=1}^{n}\left\langle dL\circ
S,\frac{\partial}{\partial v^{i}}\right\rangle\Psi^{i},$

and as $S({\Gamma})=0$ and $S(\partial/\partial v^{i})=0$, $\forall
i=1,\ldots,n$, we find the coordinate expression of $\Theta_{L}$,

$\Theta_{L}=L\,dt+\sum_{i=1}^{n}\frac{\partial L}{\partial
v^{i}}\,\theta^{i}.$ (5.9)

The exterior differential $d\Theta_{L}$ is used to define the 2-form
$\Omega_{L}=-d\Theta_{L}$, called Cartan 2-form. When $\Omega_{L}$ is of rank
$2n$ the Lagrange function is called regular. In this case, there is only one
vector field $\Gamma$ such that $i(\Gamma)\Omega_{L}=0$ and $\langle
dt,\Gamma\rangle=1$, which is the dynamical vector field defined by the
regular function $L$.

The coordinate expression of the Cartan 2-form is:

$\Omega_{L}=-d\Theta_{L}=\sum_{i,j=1}^{n}\left(-A_{ij}\theta^{i}\wedge\theta^{j}+W_{ij}\theta^{i}\wedge\Psi^{j}\right),$
(5.10)

with the functions $A_{ij}$ and $W_{ij}$ being given by

$\qquad A_{ij}=\frac{\partial^{2}L}{\partial x^{i}\partial v^{j}},\quad
W_{ij}=\frac{\partial^{2}L}{\partial v^{i}\partial v^{j}},\quad
i,j=1,\ldots,n.$ (5.11)

We recall that the dynamical vector field $\Gamma$ satisfies the two
conditions $i(\Gamma)\theta^{i}=0$ and $i(\Gamma)\Psi^{j}=0$ and, because of
this, it satisfies also the equation $i(\Gamma)\Omega_{L}=0$. Furthermore,
$(\Omega_{L})^{\wedge n}$ is a real multiple of $(\det
W_{ij})\,\theta^{1}\wedge\ldots\wedge\theta^{n}\wedge\Psi^{1}\wedge\ldots\wedge\Psi^{n}$.
Note also that as $dt\wedge\theta^{i}=dt\wedge dx^{i}$ and
$dt\wedge\Psi^{i}=dt\wedge dv^{i}$, we have that $dt\wedge(\Omega_{L})^{\wedge
n}=(\det W_{ij})\bar{\Omega}$ where $\bar{\Omega}$ denotes the volume form in
the evolution space $\bar{\Omega}=dt\wedge dx^{1}\wedge\cdots\wedge
dx^{n}\wedge dv^{1}\wedge\cdots\wedge dv^{n}$. But as
$\mathcal{L}_{\Gamma}\Omega_{L}=d(i(\Gamma)\Omega_{L})=0$ and
$\mathcal{L}_{\Gamma}dt=0$, we have that the vector field $\Gamma$ is
divergence-free with respect to the volume form $dt\wedge(\Omega_{L})^{\wedge
n}$, Consequently the determinant of the matrix with elements $W_{ij}$ is a
Jacobi multiplier with respect to the volume form in the evolution space
$\bar{\Omega}=dt\wedge dx^{1}\wedge\cdots\wedge dx^{n}\wedge
dv^{1}\wedge\cdots\wedge dv^{n}$.

Therefore, this allows us to find as a particular case the result of [14],
because if the vector field with coordinate expression

$Y=Y^{0}(t,x,\dot{x})\frac{\partial}{\partial
t}+\sum_{i=1}^{n}\left(Y^{i}(t,x,\dot{x})\frac{\partial}{\partial
x^{i}}+\bar{Y}^{i}(t,x,\dot{x})\frac{\partial}{\partial v^{i}}\right),$

is a symmetry of the vector field $\Gamma$ derived from a Lagrangian $L$,
then, according to (2.5) the function

$I={\rm div\,}Y+Y(\log W)$

is such that $\Gamma(I)=0$, where in this case the divergence should be
computed with respect to the volume form $\Omega$:

${\rm div\,}Y=\frac{\partial Y^{0}}{\partial
t}+\sum_{i=1}^{n}\left(\frac{\partial Y^{i}}{\partial
x^{i}}+\frac{\partial\bar{Y}^{i}}{\partial v^{i}}\right).$

Having in mind that

$[Y,\Gamma]=-\Gamma(Y^{0})\frac{\partial}{\partial
t}+\sum_{i=1}^{n}\left((\bar{Y}^{i}-\Gamma(Y^{i}))\frac{\partial}{\partial
x^{i}}+(Y(F^{i})-\Gamma(\bar{Y}^{i}))\frac{\partial}{\partial v^{i}}\right),$

we see that $[Y,\Gamma]=0$ if and only if $\Gamma(Y^{0})=0$,
$\bar{Y}^{i}=\Gamma(Y^{i})$ and $\Gamma(\Gamma(Y^{i}))=Y(F^{i})$. Furthermore,
$Y$ is an infinitesimal symmetry of the distribution generated by $\Gamma$,
i.e. there exists a function $h\in C^{\infty}(M)$ such that
$[Y,\Gamma]=h\,\Gamma$, if and only if $h=-\Gamma(Y^{0})$,
$\bar{Y}^{i}=\Gamma(Y^{i})+h\,v^{i}$ and
$Y(F^{i})=\Gamma(\bar{Y}^{i})+h\,F^{i}$. Therefore we obtain

$\bar{Y}^{i}=\Gamma(Y^{i})-\Gamma(Y^{0})\,v^{i},$

and

$Y(F^{i})=\Gamma(\Gamma(Y^{i}))-2\Gamma(Y^{0})\,F^{i}-\Gamma(\Gamma(Y^{0}))v^{i}.$

The two particular cases $Y^{0}\equiv 0$ and $Y^{i}\equiv 0$ were studied in
[22].

Particularly interesting cases are those for which the functions $Y^{0}$ and
$Y^{i}$ do not depend on velocities because in this case $\Gamma(Y^{0})$ and
$\Gamma(Y^{i})$ do not depend on the choice of the NSODE vector field $\Gamma$
but only of its NSODE character. On the other side, these vector fields play a
relevant role in the geometric formulation of Noether theorem, where the
symmetry of the dynamical vector field is analysed in terms of symmetry under
this kind of vector fields of the Lagrangian itself. We first recall some
basic notions of geometry of the evolution space.

In the study of nonautonomous systems of second-order differential equations
the contact 1-forms $\theta^{i}$ play a relevant role, because they can be
used to characterize curves $\gamma$ in $\mathbb{R}\times TQ\to\mathbb{R}$
that are prolongations of sections of the trivial bundle $\mathbb{R}\times Q$
over $\mathbb{R}$: In fact, such curves are characterized by
$\gamma^{*}\theta^{i}=0$, for $i=1,\ldots,n$.

Accordingly, given a vector field $X\in\mathfrak{X}(\mathbb{R}\times Q)$ there
exists one vector field $X^{(1)}\in\mathfrak{X}(\mathbb{R}\times TQ)$ such
that it is projectable on $X$ and preserves the contact distribution. If the
coordinate expression of $X\in\mathfrak{X}(\mathbb{R}\times Q)$ is

$X=X^{0}(t,x)\frac{\partial}{\partial
t}+\sum_{i=1}^{n}\left(X^{i}(t,x)\frac{\partial}{\partial x^{i}}\right),$
(5.12)

then a vector field $X^{(1)}\in{\hbox{{\frak X}}}(\mathbb{R}\times TQ)$
projectable on $X$, must be of the form

$X^{(1)}=X^{0}(t,x)\frac{\partial}{\partial
t}+\sum_{i=1}^{n}\left(X^{i}(t,x)\frac{\partial}{\partial
x^{i}}+\bar{X}^{i}(t,x,v)\frac{\partial}{\partial v^{i}}\right),$ (5.13)

where, choosing an arbitrary NSODE vector field and the corresponding
associated bases of 1-forms and vector fields, as

$\mathcal{L}_{X^{(1)}}\theta^{i}=\sum_{j=1}^{n}\left(\frac{\partial
X^{i}}{\partial x^{j}}-v^{i}\frac{\partial X^{0}}{\partial
x^{j}}\right)\,\theta^{j}+(\Gamma X^{i}-\bar{X}^{i}-v^{i}\,\Gamma X^{0})\,dt,$

this shows that in order to the set of contact 1-forms be preserved by
$X^{(1)}$, it must be of the form

$\bar{X}^{i}=\Gamma(X^{i})-v^{i}\,\Gamma(X^{0}),\quad i=1,\ldots,n.$ (5.14)

Note that the last relations are valid for any NSODE $\Gamma$.

Next we consider a particular case that contains some of the results obtained
in [22].

In the particular case of a vector field $X^{(1)}$ that is a symmetry of the
one-dimensional distribution generated by a NSODE vector field $\Gamma$
determined by a Lagrangian $L$, as $\det W$ is a Jacobi multiplier for
$\Gamma$ the Hojman constant of the motion according to (2.6) is given by

$I=2\sum_{i=1}^{n}\left(\frac{\partial X^{i}}{\partial
x^{i}}-v^{i}\frac{\partial X^{0}}{\partial
x^{i}}\right)-n\,\Gamma(X^{0})+X^{(1)}(\log(\det W)),$ (5.15)

because

${\rm div\,}(X^{(1)})=\frac{\partial X^{0}}{\partial
t}+\sum_{i=1}^{n}\left(\frac{\partial X^{i}}{\partial
x^{i}}+\frac{\partial\bar{X}^{i}}{\partial v^{i}}\right),$

and as $X^{0}$ and $X^{i}$ do not depend on the velocities,

$\frac{\partial(\Gamma(X^{0}))}{\partial v^{i}}=\frac{\partial X^{0}}{\partial
x^{i}},\qquad\frac{\partial(\Gamma(X^{i}))}{\partial v^{i}}=\frac{\partial
X^{i}}{\partial x^{i}},\,$

and then,

${\rm div\,}(X^{(1)})=\frac{\partial X^{0}}{\partial
t}+\sum_{i=1}^{n}\left(\frac{\partial X^{i}}{\partial
x^{i}}+\frac{\partial(\Gamma(X^{i})-v^{i}\,\Gamma(X^{0}))}{\partial
v^{i}}\right),$

and therefore

${\rm div\,}(X^{(1)})=\frac{\partial X^{0}}{\partial
t}+\sum_{i=1}^{n}\left(2\frac{\partial X^{i}}{\partial
x^{i}}-v^{i}\frac{\partial X^{0}}{\partial x^{i}}-\Gamma(X^{0})\right),$

which can be rewritten as

${\rm div\,}(X^{(1)})=2\sum_{i=1}^{n}\left(\frac{\partial X^{i}}{\partial
x^{i}}-v^{i}\frac{\partial X^{0}}{\partial x^{i}}\right)-(n-1)\Gamma(X^{0}),$

and, as $\det W$ is a Jacobi multiplier, the constant of the motion
corresponding to (4.3) is (5.15), in agreement with the above mentioned
results in [22].

##  6 Conclusions and outlook

It has been shown how the use of basic results of the differential geometry
concerning the concept of divergence of a vector field, as well as some
properties relating Lie derivatives of forms with respect to couples of vector
fields, allows us to establish in Theorem 1 a method to find a first-integral
for a divergence-free vector field when a symmetry vector field is known. The
theory is generalised to include a more general vector field when a Jacobi
multiplier for such given vector field is known in Theorem 2. Theorem 3
presents the extension of the results of these Theorems when a normalizer
vector field instead of a symmetry vector is known. This generalization is
very relevant for the case of a nonautonomous system. The main interest of
these Theorems is that they present in an intrinsic way a general result that
contains as particular cases many expressions of first-integrals or constants
of motion derived by ad-hoc methods in many branches of physics. This makes
very interesting this unifying way of finding them using geometric tools of
modern differential geometry. Several examples of applications to general
Hamiltonian systems, second order differential equations and Lagrangian
systems, both in the autonomous and the nonautonomous cases, have also been
exhibited.

## References

  * [1] D.G. Currie and E.J. Saletan, $q$-equivalent particle Hamiltonians. The classical one-dimensional case, J. Math. Phys. 7, 967–974 (1966).
  * [2] S. Hojman and H. Harleston, Equivalent Lagrangians: multidimensional case, J. Math. Phys. 22, 1414–1419 (1981).
  * [3] J.F. Cariñena and L.A. Ibort, Non-Noether constants of motion, J. Phys. A: Math. Gen. 16, 1–7 (1983).
  * [4] S. Hojman, A new conservation law constructed without using either Lagrangians or Hamiltonians, J. Phys. A. Math. Gen. 25, L291–L295 (1992).
  * [5] F. González-Gascón, Geometric foundations of a new conservation law discovered by Hojman, J. Phys. A: Math. Gen. 27, L 59–L 60 (1994).
  * [6] F. Darabi, M. Golmohammadi, and A. Rezaei-Aghdam, FRW string cosmological solutions via Hojman symmetry, Int. J. Geom. Methods Mod. Phys. 17, 2050175 (2020).
  * [7] F. Darabi, M. Golmohammadi, and A. Rezaei-Aghdam, Generalized (2 + 1)-dimensional BTZ black holes via Hojman symmetry, arXiv:2010.08424v2 (2021).
  * [8] H. Wei, Y-N. Zhou, H-Y. Lii and X-B. Zou, Hojman symmetry in $f(T)$ theory, Astrophys. Space Sci. 360, 6 (2015).
  * [9] H. Wei, Y-N. Zhou, H-Y. Lii and X-B. Zou, Exact cosmological solutions of $f(R)$ theories via Hojman symmetry, Nucl. Phys. B 903, 132–149 (2016).
  * [10] S. Capozziello and M. Roshanc, Exact cosmological solutions from Hojman conservation quantities, Phy. Lett. B 726, 471–480 (2013).
  * [11] M.C. Paolella and S. Capozziello, Hojman symmetry approach for scalar-tensor cosmology, Phy. Lett. A 379, 1304–1308 (2015).
  * [12] A. Paliathanasis, P.G.L. Leach and S. Capozziello, On the Hojman conservation quantities in Cosmology, Phy. Lett. B 755, 8–12 (2016).
  * [13] F. González-Gascón, Notes on the Connection between the Symmetries and the First Integrals of Dynamical Systems, Lett. Nuovo Cim. 19, 366–368 (1977).
  * [14] M. Lutzky, Remarks on a recent theorem about conserved quantities, J. Phys. A: Math. Gen. 28, L 637–L 638 (1995).
  * [15] J.F. Cariñena and P. Santos, Jacobi Multipliers and Hamel’s formalism, J. Phys. A: Math. Theor. 54, 225203 (2021).
  * [16] J.F. Cariñena, A. Ibort, G. Marmo and G. Morandi, Geometry from Dynamics: Classical and Quantum. Springer, Dordrecht, 2015.
  * [17] F. Mei, Non-Noether conserved quantity for differential equations of motion in the phase space, Chin, Sci. Bull. 47, 2049–2050 (2002).
  * [18] S-L. Gu and K-X, Wei, Study on The Symmetry and Conserved Quantities for Hamilton Systems, International Conference on Logistics Engineering, Management and Computer Science (LEMCS 2014), Advances in Intelligent Systems Research Series, p. 798–801 (2014).
  * [19] M. Crampin, Constants of the Motion in Lagrangian Mechanics, Int. J. Theor. Phys. 16, 741–754 (1977).
  * [20] J.F. Cariñena, J. Fernández–Núñez and E. Martínez, A geometric approach to Noether’s Second Theorem in time-dependent Lagrangian Mechanics Lett. Math. Phys. 23, 51–63 (1991).
  * [21] J.F. Cariñena, J. Fernández–Núñez and E. Martínez, Noether’s theorem in time-dependent Lagrangian Mechanics, Rep. Math. Phys. 31, 189–203 (1992).
  * [22] H-B. Zhang and L-Q. Chen, The Unified Form of Hojman’s Conservation Law and Lutzky’s Conservation Law, J. Phys. Soc. Japan 74, 905–909 (2005).
  * [23] F. González-Gascón, Divergence-free vector fields and integration via quadratures, Phys. Lett. A 225, 269–273 (1996).
  * [24] M. Lutzky, Conserved quantities and velocity-dependent symmetries in Lagrangian dynamics, Int. J. Nonlinear Mechanics 33, 393–396 (1998).
  * [25] W.-L. Huang and J.-L. Cai, Inverse Problems of Mei Symmetry for Nonholonomic Systems with Variable Mass, J. Mech. 31, 515–523 (2015).
  * [26] S. Hojman, Symmetries of Lagrangians and of their equations of motion, J. Phys. A, Math. Gen. 17, 2399-2412. (1984).
  * [27] S. Hojman, Interpretation of symmetry transformations, J. Phys. A, Math. Gen. 17, L521–L525 (1984).
  * [28] S-K. Luo, New Types of the Lie Symmetries and Conserved Quantities for a Relativistic Hamiltonian System, Chinese Phys. Lett. 20, 597–599 (2003).
  * [29] Y.-L. Han, X.-X. Wang, M.-L. Zhang and L.-Q. Jia, Special Lie symmetry and Hojman conserved quantity of Appell equations for a Chetaev nonholonomic system, Nonlinear Dyn. 73, 357–361 (2013).
  * [30] Y.-L. Han, X.-X. Wang, M.-L. Zhang and L.-Q. Jia, Lie Symmetry and Approximate Hojman Conserved Quantity of Lagrange Equatio ns for a Weakly Nonholonomic System, J. Mech. 30, 21–27 (2014).
  * [31] J-L. Fu, L-Q. Chen, S. Jiménez and Y-F. Tang, Non-Noether symmetries and Lutzky conserved quantities for mechanico-electrical systems, Phys. Lett. A 358, 5–10 (2006).