# Spin dynamics of the quantum dipolar magnet Yb3Ga5O12 in an external field

E. Lhotel  Institut Néel, CNRS & Univ. Grenoble Alpes, 38000 Grenoble, France
L. Mangin-Thro  Institut Laue Langevin, 38042 Grenoble, France  E. Ressouche
Univ. Grenoble Alpes, CEA, IRIG, MEM, MDN, 38000 Grenoble, France  P. Steffens
Institut Laue Langevin, 38042 Grenoble, France  E. Bichaud  Univ. Grenoble
Alpes, Grenoble INP, CEA, IRIG, PHELIQS, 38000 Grenoble, France  G. Knebel
Univ. Grenoble Alpes, Grenoble INP, CEA, IRIG, PHELIQS, 38000 Grenoble, France
J.-P. Brison  Univ. Grenoble Alpes, Grenoble INP, CEA, IRIG, PHELIQS, 38000
Grenoble, France  C. Marin  Univ. Grenoble Alpes, Grenoble INP, CEA, IRIG,
PHELIQS, 38000 Grenoble, France  S. Raymond<EMAIL_ADDRESS>Univ. Grenoble
Alpes, CEA, IRIG, MEM, MDN, 38000 Grenoble, France  M. E. Zhitomirsky
<EMAIL_ADDRESS>Univ. Grenoble Alpes, Grenoble INP, CEA, IRIG,
PHELIQS, 38000 Grenoble, France

###### Abstract

We investigate ytterbium gallium garnet Yb3Ga5O12 in the paramagnetic phase
above the supposed magnetic transition at $T_{\lambda}\approx 54$ mK. Our
study combines susceptibility and specific heat measurements with neutron
scattering experiments and theoretical calculations. Below 500 mK, the elastic
neutron response is strongly peaked in the momentum space. Along with that the
inelastic spectrum develops flat excitation modes. In magnetic field, the
lowest energy branch follows a Zeeman shift in accordance with the field-
dependent specific heat data. An intermediate state with spin canting away
from the field direction is evidenced in small magnetic fields. In the field
of 2 T, the total magnetization almost saturates and the measured excitation
spectrum is well reproduced by the spin-wave calculations taking into account
solely the dipole-dipole interactions. The small positive Curie-Weiss
temperature derived from the susceptibility measurements is also accounted for
by the dipole spin model. Altogether, our results suggest that Yb3Ga5O12 is a
quantum dipolar magnet.

##  I Introduction

Geometrical frustration typically suppresses conventional magnetic ordering
and may stabilize exotic phases with unconventional spin excitations [1, 2].
Apart from fundamental interest, the delayed magnetic ordering in conjunction
with a large unfrozen entropy makes frustrated materials to be good candidates
for the low-temperature magnetic cooling [3]. The developing space
applications and increasing costs of helium fuel a continuing search for new
refrigerant materials for the adiabatic demagnetization refrigeration (ADR) in
the 0.1–4 K temperature range, which can outperform the presently used
paramagnetic salts [4]. In the past few years there have been experimental
reports of the enhanced magnetocaloric effect observed for pyrochlores [5, 6,
7], intermetallics [8, 9, 10], and various Yb compounds [9, 10, 11, 12].

In this work, we investigate ytterbium gallium garnet Yb3Ga5O12, which
exhibits an enhanced magnetocaloric effect below 2 K [11]. There is also a
reborn interest in the rare-earth gallium garnets $R_{3}$Ga5O12 in general.
The recent experimental studies focus on $R=\rm Er$ [13], $R=\rm Tb$ [14, 15],
and $R=\rm Dy$ [16] compounds, whereas gadolinium gallium garnet was
extensively studied in the past two decades, see, for example, [17] and
references therein.

In the cubic garnet crystals, space group $Ia\bar{3}d$, the rare-earth ions
form two inter-penetrating lattices of corner-sharing triangles often called
hyperkagome lattices. The hyperkagome lattice strongly frustrates the nearest-
neighbor antiferromagnetic exchange keeping a large degree of frustration in
the presence of anisotropy and dipolar interactions as well. Positions of the
crystal field levels of Yb3+ ions were measured by Buchanan et al. for several
ytterbium garnets using infrared spectroscopy [18]. The lowest $J=7/2$
multiplet is split into four Kramers doublets by the local orthorhombic
crystal field. In Yb3Ga5O12, the distance between the lowest doublet and the
next excited level is about 67 meV. Hence, at temperatures below 10 K, one can
safely neglect the excited levels and interpret magnetic properties in terms
of the effective $S=1/2$ degrees of freedom corresponding to the ground state
doublet. The magnetic entropy variation $\Delta S_{m}=R\ln 2$ determined in
the specific heat measurements between 44 mK and 4 K [19] is fully consistent
with this conclusion. The $g$-tensor of the ground-state doublet is only
weakly anisotropic as follows from EPR measurements [20] and crystal-field
calculations [21].

Below 4 K, the magnetic susceptibility of Yb3Ga5O12 follows the Curie-Weiss
law with a small Curie-Weiss temperature $\theta_{p}=45$ mK [19]. In addition,
the specific heat exhibits a lambda-peak anomaly at $T_{\lambda}=54$ mK
preceded by a broad maximum at 200 mK. These features were tentatively
associated with a long-range order and short-range spin correlations,
respectively [19]. However, Mössbauer spectroscopy measurements did not detect
any hyperfine splitting associated with a possible magnetic ordering below
$T_{\lambda}$ [22]. $\mu$SR experiments did not observe signs of a long-range
order as well pointing instead to the reinforcement of the dynamical
correlations below $T_{\lambda}$ [23]. The diffuse magnetic scattering in
Yb3Ga5O12 has been recently measured in neutron experiments [24] and
interpreted in terms of the exotic director state associated with ten-site
loops by analogy with a similar suggestion for Gd3Ga5O12 [17]. The
experimental studies in Ref. [24] have been complemented by simulations of the
exchange-dipolar Ising model. Thus, despite several decades of research, the
basic magnetic properties of ytterbium gallium garnet remain poorly understood
both in terms of possible magnetic states [25] as well as underlying
microscopic interactions. In the present paper, we combine results of bulk
measurements of single crystals of Yb3Ga5O12 with neutron scattering
experiments in an applied magnetic field and theoretical calculations to get a
further insight into the microscopic properties of this interesting frustrated
material. The combination of complementary techniques allows us to
quantitatively conclude on a dominant role of the dipolar interactions both
for the static and dynamic properties of this interesting spin-1/2 hyperkagome
material.

##  II Experimental details

###  II.1 Sample preparation

Single crystals of Yb3Ga5O12 were grown by the floating zone technique using a
commercial optical furnace (Crystal System Inc.) equipped with 1 kW halogen
lamps. The feed rods were prepared from high-purity oxides Yb2O3 (99.99%) and
Ga2O3 (99.99%) mixed and heated up to 1150∘C. After sintering, the rods were
pressed isostatically at 1500 bar and sintered again at 1400∘C in air. Crystal
growth conditions were optimized under argon pressure (4.5 bar) with a gas
flow of 0.5 L/min at the growth rate of 2 mm/h. The feed and crystallized rods
have been rotated in opposite directions at the rate of 30 rounds per minute.
Crystals of Yb3Ga5O12 prepared with this procedure grow in a direction very
close to a four-fold axis of the cubic $Ia{\bar{3}}d$ phase. We observe that
rods are white transparent when grown under air or oxygen (oxidizing)
conditions or slightly pink under argon (reducing) atmosphere. No foreign
phase was detected in the $X$-ray powder diffraction experiments. The checks
of crystal quality and the crystallographic orientation were done using the
Laue diffraction technique. It should be mentioned that the crystal color can
be reversibly changed from white to pink according to the post growth heat
treatment in a specific atmosphere. No significant modifications of the basic
physical properties (magnetization, specific heat) were observed between both
states.

Growth of single crystals was complemented by preparation of a high-density
ceramic of Yb3Ga5O12 suitable for a direct use in an ADR prototype [11].
Ceramic samples were made from the same starting powders of Yb2O3 and Ga2O3.
An extensive investigation was carried out using a planetary crusher to
perform fine grinding followed by laser granulometry to measure a powder-
diameter distribution, and high-temperature dilatometry to determine optimal
temperatures for synthesis and sintering. After all, the starting materials
were first treated at 1050∘C before being pressed isostatically under 2500 bar
and, then, sintered at 1700∘C. The ceramic samples were checked by $X$-ray
powder diffraction, showing the same cubic $Ia\bar{3}d$ phase. No traces of
foreign phases have been detected. Following this optimized route, the density
of ceramic samples reaches 96(1)% of the bulk crystal density showing no open
porosity as evidenced by scanning electronic microscopy.

###  II.2 Methods

The specific heat was measured on a commercial Quantum Design Physical
Property Measurement System with the 3He insert using the heat-pulse semi-
adiabatic technique. Two types of samples were investigated: a 0.87 mg ceramic
sample in fields up to 8 T and a 17.25 mg single crystal (white color) with
the field applied along the [001] direction in fields up to 3 T.

Magnetization was measured using the home-made low temperature SQUID
magnetometer developed at the Institut Néel, which operates in a dilution
fridge [26]. A single crystal of Yb3Ga5O12 shaped as a parallelepiped of size
$0.93\times 0.95\times 4.37$ mm3 has been glued on a copper sample holder with
GE varnish to ensure a good thermalization. Magnetic field was applied along
the long side of the sample, which corresponds to the [001] direction. The
corresponding demagnetization factor was estimated to 1.19 (cgs units) [27].
All the data shown below are corrected from demagnetizing effects.

The neutron scattering experiments were performed at the high-flux reactor of
the Institut Laue-Langevin, Grenoble. Diffraction measurements were done on
the two-axis thermal-neutron diffractometer D23 equipped with a lifting
detector. A copper monochromator provides an unpolarized beam with the
wavelength of 1.276 Å. The diffuse scattering was measured on the D7 beam line
for the polarization along the vertical $z$-axis and the wavelength of the
incident beam being 4.8 Å. A four-to-one time ratio was used between the spin-
flip and non spin-flip channel measurements. The data were collected in a 120∘
range, with steps of 1∘ for two positions of the detector banks, resulting in
0.5∘ between each measured position. Standard vanadium, quartz and empty
corrections have been performed for taking into account the detector
efficiency, the polarization efficiency and the background scattering [28].
The flipping ratio measured in the direct beam was in the range 20-29 during
all the experiment. However the flipping ratio was larger at higher
temperatures (20-27 at 50 mK and 25-29 at 5 K).

The inelastic neutron-scattering (INS) experiments were carried out on the
three-axis spectrometer ThALES. The instrument was used in the classical
single detector mode and set up in the W configuration. The incident beam was
provided by the Si monochromator with a velocity selector in place. A Be
filter was put in the scattered beam and a radial collimator was placed
between the graphite analyzer and the detector. Most of the data were taken at
fixed $k_{f}=1.2$ Å-1 (corresponding to 5.2 Å).

The white-color crystal of Yb3Ga5O12 used in the neutron experiments has a
cylindrical shape (diameter 8 mm, length 27 mm) with the long axis parallel to
the [001] direction. The sample was embedded in a copper foil glued with
hydrogen-free fluorinated Fomblin oil, tightened with Teflon tape to ensure
good thermalization, and mounted in a dilution insert. For D23, a smaller part
cut from this rod was used. The scattering plane was selected to be (1 0 0)-(0
1 0). For D23 and ThALES, the sample was inserted inside 6 T and 2.5 T
vertical magnets respectively, so that the field was applied along the [001]
direction. In our experiments, the lowest temperature indicated on the
dilution fridge thermometer was 50 mK. Still, due to the distance between the
sample and the dilution mixing chamber as well as to a heating effect from the
neutron beam, we assume that our measurements were performed above 54 mK,
which is the expected transition temperature for Yb3Ga5O12 [19]. For the sake
of simplicity, we nevertheless use the thermometer temperatures in the
description of our results. In the following, the scattering vectors are
expressed as $\textbf{Q}=(Q_{h},Q_{k},Q_{l})$ and reciprocal space directions
are labelled by $(h\ k\ l)$.

###  II.3 Magnetic characterization

The static magnetic susceptibility was measured in a weak applied field
$\mu_{0}H_{e}=0.01$ T to insure a linear response in the whole range of
temperatures down to 70 mK. The susceptibility was calculated by dividing a
uniform magnetization to an internal field $H_{i}$, which was obtained from
$H_{e}$ by applying the demagnetization correction at each temperature. The
inverse susceptibility can be fitted in the 1–4.2 K range to a Curie-Weiss law
$\chi^{-1}=(T-\theta_{p})/C$, where $\theta_{p}$ is the Curie-Weiss
temperature and $C$ is the Curie constant (see Fig. 1). The fitting yields
$\theta_{p}=97$ mK and $C=1.12$ emu/mol Yb. The Curie constant corresponds to
an effective moment of $3~{}\mu_{\rm B}$ for an Yb3+ ion in agreement with the
early measurements by Fillipi et al. [19]. Our value for the Curie-Weiss
temperature is somewhat larger than their estimate $\theta_{p}\simeq 45$ mK,
though both agree on its ferromagnetic sign. In contrast, a negative
(antiferromagnetic) $\theta_{p}$ has been recently reported in Ref. 24. Note,
that accounting for the Van Vleck contribution, as was done in Ref. 19, only
slightly affects the values of $C$ and $\theta_{p}$.

Figure 1: Inverse magnetic susceptibility $\chi^{-1}=H_{i}/M$ as a function of
temperature measured in the applied field $\mu_{0}H_{e}=0.01$ T oriented
parallel to the [001] axis. The solid line is a fit to the Curie-Weiss law.

The saturated magnetization measured at 120 mK and 8 T is found to be 1.76
$\mu_{\rm B}$, corresponding to an effective $g=3.52$, in agreement with the
previous reports and the crystal field calculations [19, 20, 21].

##  III Diffuse scattering

The magnetic scattering measured on D7 in the spin flip channel with
polarization along the $c$-axis is shown in Fig. 2 for four different
temperatures 50 mK, 200 mK, 500 mK and 5 K. The measured magnetic signal
corresponds to the scattering obtained for an energy integration of the
scattering function in the range [-$k_{B}T$, $E_{i}$], with $E_{i}$=3.5 meV.

Magnetic correlations are observed along the (1 1 0) direction and are peaked
at $Q\approx 0.36$ Å-1 (feature A) and $Q\approx 1.1$ Å-1 (features B). They
are washed out progressively with temperature and a $\bf{Q}$-independent
weaker signal is observed at 5 K. In addition a small signal with the same
characteristics in temperature is observed for $Q\approx 1.95$ Å-1 (in the
directions (3 7 0) and (7 3 0); features C). The observation of Bragg peaks at
high momentum transfer is ascribed to the inaccuracy of the polarization
correction which happens for these nuclear Bragg peaks that have a tremendous
intensity in the non spin flip channel.

Figure 2: Color-coded intensity plot of the magnetic scattering measured on D7
at (a) 50 mK, (b) 200 mK, (c) 500 mK, and (d) 5 K. The labels A, B and C in
panel (a) mark the diffuse scattering features described in the text and the
stars indicate the location of the momentum transfer ${\bf Q}_{1}=(3/2,3/2,0)$
and ${\bf Q}_{2}=(3,0,0)$ investigated by INS. The dashed lines indicate the
longitudinal and transverse cuts shown in Fig. 3.

Figure 3(a) shows the peak obtained by performing a cut of the feature A in
the direction (1 -1 0) for the four measured temperatures. The peaks are
fitted by Lorentzian lineshapes and the resulting temperature dependence of
the amplitude $A$, and Half Width at Half Maximum (HWHM) $\kappa$, are shown
in Fig. 3(b). The signal clearly rises up upon cooling, and the correlation
length increases below 200 mK. At 50 mK, the obtained HWHM is 0.12 r.l.u.
(along (1 -1 0)), which corresponds to a correlation length, $\xi\propto
1/\kappa$, of about 11 Å, i.e. to a distance between several Yb atoms in the
cell (the minimum Yb-Yb distance being 3.74 Å).

Figure 3: (a) Momentum cut across feature A along the (1 -1 0) direction for
several temperatures, lines show the Lorentzian fits. (b) Temperature
variation of the amplitude, $A$, and HWHM, $\kappa$, of the peaks shown in
(a). Lines are guides to the eye. (c) Constant energy spectra performed on
ThALES along ${\bf Q}=(Q_{h}$, $Q_{h}$, 0) at $E=0$ meV in zero field at
$T=65$ mK and 1.4 K, subtracted from the 65 mK data for $\mu_{0}H=2$ T. The
same data without subtraction are shown in the inset (in the 2 T data, the
missing points in the figure correspond to the (1, 1, 0) magnetic Bragg peak
induced under field). Dashed lines are guides to the eyes. In the main panel,
the full line corresponds to the same longitudinal cut obtained from the D7
data at 50 mK.

The elastic scattering was measured on ThALES with an unpolarized neutron
setup. Here, the measurement is made for zero energy transfer, $E=0$, and the
integration in energy is performed over the resolution function, which Full
Width at Half Maximum (FWHM) equals 0.039 meV as measured on the incoherent
signal. The obtained scans along ${\bf Q}=(Q_{h},\ Q_{h},\ 0)$ at 65 mK and
1.4 K are shown in Fig. 3(c).

A modulation is observed at 65 mK, which is suppressed at 1.4 K. The observed
signal at $E=0$ on ThALES corresponds to the features A and B highlighted on
the D7 data as demonstrated by the similar peak positions in the figure. This
indicates that the diffuse response is essentially dominated by true elastic
or very low energy processes within the neutron time-scale. It is to be noted
that the sample base temperature reached on ThALES is undoubtly higher than
the one reached on D7 given the very high flux of cold neutron on the former
instrument and this could explain the difference of lineshapes obtained owing
to the strong temperature dependence of the signal below 200 mK.

##  IV Zero field excitation spectrum

The spin dynamics has been investigated along the high-symmetry lines and
points of the Brillouin zone. Emphasis was put on specific wave-vectors at
which the magnetic diffuse scattering described in the previous section was
found to have its maximum or minimum intensity.

Representative low energy spectra are shown in Fig. 4 for ${\bf
Q}=(3/2,3/2,0)$, panels (a) and (b), and ${\bf Q}=(3,0,0)$, panels (c) and
(d), at $T=65$ mK and 1.4 K. Positions of these two vectors are indicated by
${\bf Q}_{1}$ and ${\bf Q}_{2}$ in Fig. 2. At the lowest temperature, aside
from the elastic line, the spectra are composed of several inelastic modes: a
sharp mode at $E_{0}$, around 0.05 meV, a broader one at $E_{1}$, around 0.1
meV whose maximum intensity is barely resolved [panels (a) and (c)] and a very
broad one at $E_{2}$, around 0.7 meV [panels (b) and (d)]. The constant
Q-scans are phenomenologically described at low energy, panels (a) and (b), by
a Gaussian lineshape, where both Stokes and anti-Stokes peaks together with
the temperature factor are taken into account. A constant background and an
elastic line with a fixed profile are also included (See Appendix A).

Figure 4: Energy spectra at $T=65$ mK and 1.4 K in zero field: (a) and (b) for
the scattering vector ${\bf Q}=(3/2,\ 3/2,\ 0)$ and (c) and (d) for ${\bf
Q}=(3,\ 0,\ 0)$. Lines are fits as explained in the text.

The fits shown in Fig. 4 are obtained with $E_{0}=0.045(6)$ meV and
$E_{1}=0.100(6)$ meV and fixed FWHM of 0.059 meV and 0.094 meV respectively,
which are about 2 and 3 times the incoherent width. The higher energy part is
conveniently described by a broad Gaussian peak centered at $E_{2}=0.70(5)$
meV with a FWHM of 0.190(5) meV on a slopping background (See panels (b) and
(d) of Fig. 4). At 1.4 K, the overall intensity is reduced in both the elastic
and inelastic channels and a broad quasielastic-like signal replaces the well
defined excitations. It is thus impossible to describe in a unique manner the
lowest energy modes. The fits shown in Fig. 4 are indicative only and are
obtained keeping the peak positions $E_{0}$ and $E_{1}$ to their low
temperature values, while the peak widths and intensities are free. On the
other hand, the spectral weight of the $E_{2}$ mode clearly moves to lower
energies.

Figure 5: (a) Constant energy scans performed along ${\bf Q}=(Q_{h},\ Q_{h},\
0)$ for the energy transfer $E=0.05$ meV at $T=65$ mK and 1.4 K. (b) Energy
spectra at ${\bf Q}=(Q_{h},\ Q_{h},\ 0)$ for $Q_{h}=1/2$, 1 and $3/2$ and at
${\bf Q}=(0.54,\ 2.05,\ 0)$ at 65 mK. Lines are fits performed using the same
procedure as described for Fig. 4.

The inelastic modes are essentially dispersionless and a representative
constant energy scan performed along ${\bf Q}=(Q_{h},\ Q_{h},\ 0)$ for an
energy transfer of 0.05 meV is shown in the panel (a) of Fig. 5. A broad
intensity modulation with Q is observed at low temperature, with a kink at
$Q_{h}\approx 0.75$ r.l.u corresponding to a scattering vector magnitude of
0.55 Å-1 (arrow in Fig. 5(a)). The dispersionless nature of the excitations is
further confirmed by several constant Q scans performed along this line for
$Q_{h}=1/2$, 1 and $3/2$ shown in Fig. 5(b) together with a scan performed at
an arbitrary position $\textbf{Q}=(0.54,\ 2.05,\ 0)$ having the same Q modulus
as $\textbf{Q}=(3/2,\ 3/2,\ 0)$. Here again, only an intensity modulation is
present, while no clear dispersion in energy could be detected. In the
remaining part of the paper, we focus on the response at the scattering
vectors $\textbf{Q}=(3/2,3/2,0)$ and $\textbf{Q}=(3,0,0)$.

##  V Magnetic response

Figure 6: Energy spectra at ${\bf Q}=(3/2,3/2,0)$ at $T=65$ mK and various
magnetic fields ${\bf H}\parallel[001]$ with the intensity on a log scale (a)
and on a linear scale (b). The inset shows a magnified view on the 0–0.6 meV
energy range. Lines in panel (b) are fits described in the text.

The evolution of the excitation spectrum in an applied field is shown in Fig.
6 for $\textbf{Q}=(3/2,3/2,0)$ at $T=65$ mK. Figure 6(a) presents the
intensity on a log scale to emphasize the presence of three inelastic modes.
Their positions change with field as indicated by different arrows. The two
lowest energy modes at $E_{0}$ and $E_{1}$ move to higher energy, the first
one being better defined when escaping from the elastic line. The higher-
energy mode $E_{2}$, initially at 0.7 meV moves to lower energy and merges
with the second mode at around 1 T. Figure 6(b) shows this evolution on a
linear scale which better illustrates the behavior of the intense elastic line
and of the lowest energy inelastic mode at $E_{0}$. The inelastic line is
fully resolved at 0.5 T, but exhibits broadening with increasing field and
eventually gives rise to a double peak structure at 2 T, see the inset in Fig.
6(b). As for the elastic line, its intensity is reduced at 0.25 T compared to
zero field and does not further change up to 2 T. This observation indicates
that the diffuse magnetic scattering seen in the D7 data is fully suppressed
at 0.25 T.

Figure 7: Field dependence of physical parameters (${\bf H}\parallel[001]$):
(a) uniform magnetization; (b) intensities of the Bragg peaks $(1,1,0)$ (main
panel) and (2, 2, 0) (inset); (c) energies of the three excitation modes
obtained from the INS measurements. The full line is a linear fit, the dashed
lines are guides to the eye. Open circles show the Schottky gap obtained from
the specific heat data.

The field dependence of the peak positions is shown in Fig. 7(c). The energy
of the first mode, which is the only one well-defined in the whole field
range, is fitted with a linear variation. The intercept at 0 T is 0.032(5) meV
and the slope is 0.198(5) meV/T [29]. It is to be noted that the energy of the
first mode determined in this way is found to be lower than the one obtained
by fitting the raw data in Fig. 4, i.e. 0.032(5) versus 0.045(1) meV. The
slope conversion in terms of magnetic moments gives a value of 3.41 $\mu_{\rm
B}$. This suggests that the $E_{0}$ field dependence is associated with the
Zeeman splitting $2\mu H$, giving a moment $\mu$ of $1.70~{}\mu_{\rm B}$
consistent with the magnetization measurements below the saturation (between 2
and 8 T, the magnetization evolves from 1.70 to 1.76 $\mu_{\rm B}$).

A further evidence of the Zeeman splitting effect is found in specific heat
measurements, see Fig. 8. With increasing field, the broad peak in $C(T)$
shifts to higher temperatures. The data can be fitted to a Schottky anomaly
and the field dependence of the corresponding gap is shown in Fig. 8(b). A
linear variation can be traced up to $\mu_{0}H=4$ T. A similar field variation
of the Schottky anomaly is also found for a single crystal sample (solid
squares in the figure). The obtained gap is compared to the lowest energy
excitation in Fig. 7(c) and a good matching is obtained. Indeed, as seen
above, under magnetic field (already at 0.25 T), the diffuse magnetic elastic
scattering disappears so that the low energy correlations existing in zero
field are suppressed. It is thus natural that the specific heat Schottky
anomaly reflects the Zeeman splitting of the INS spectrum. At higher fields (6
and 8 T), the fit to a Schottky anomaly is less accurate: the specific heat
peak is narrower and of smaller intensity. This hints for a more complex
picture than a simple two level system description.

Figure 8: (a) Temperature dependence of the specific heat of a ceramic sample
for several magnetic fields. The open circles are the zero-field data by
Filippi et al. [19]. Lines are fits to the Schottky anomaly function. (b)
Field dependence of the Schottky gap $\Delta$, open circles; solid squares are
the additional single crystal data obtained for ${\bf H}\parallel[001]$. The
solid line is a linear fit to the low-field data up to 4 T.

In order to shed light on the field induced behavior of Yb3Ga5O12, neutron
diffraction measurements were conducted on D23. The lattice parameter obtained
from the orientation matrix of the single crystal sample is $a=12.19$ Å. The
focus was on the magnetic field variations of the $(2,2,0)$ and $(1,1,0)$
Bragg reflection intensities. The inset of Fig. 7(b) shows the field variation
at 65 and 500 mK of the (2, 2, 0) peak intensity. This peak, allowed by
symmetry for the 24c Yb site, probes the ferromagnetic component of the
system. Indeed, its field dependence strongly resembles the low temperature
magnetization shown in Fig. 7(a) in the low field region. At 65 mK, the signal
increases up to about 0.5 T and then slowly decreases. The field dependence of
the $(1,1,0)$ peak is shown in Fig. 7(b). (1, 1, 0) is a forbidden peak for
the 24c position and a magnetic field induced signal thus points to a breaking
of the usual uniform field induced symmetry in the achieved spin
configuration. At low temperature, the amplitude of the $(1,1,0)$ peak
increases fast and reaches a maximum at 0.15 T. This intermediate state
corresponds to a canting of some magnetic moments. The intensity decreases
above 0.15 T, but remains finite and almost stable above 1 T, which matches
the slow decrease of the (2,2,0) peak at high field. The magnetization itself
never saturates up to 8 T (not shown), a small but finite slope being observed
above 1 T. In the following, the phase at 2 T is nonetheless named saturated
for convenience.

The measured dispersion of low-lying magnetic excitations in the 2 T saturated
phase is illustrated in Fig. 9. The panels show the energy spectra for the (0
1 0) direction [30], Fig. 9(a), and the (1 1 0) direction, Fig. 9(c), measured
at $T=65$ mK. They exhibit a double peak structure with an individual peak
width of 0.05 meV (FWHM). Figures 9(b) and (d) present the combined color-
coded intensity plots for the two directions. The excitation peaks are
essentially dispersionless. The magnon band is centered near 0.4 meV and spans
an energy region of 0.1 meV in width.

Figure 9: Inelastic neutron scattering measurements at $T=60$ mK and
$\mu_{0}H=2$ T: (a) and (c) the energy spectra along ${\bf Q}=(2,\ k,\ 0)$ and
${\bf Q}=(2-h,\ 2-h,\ 0)$, respectively; (b) and (d) show the corresponding
color-coded intensity maps as functions of Q and energy transfer. In (a) and
(c) the lines are Gaussian fits. In the color plots (b) and (d), the shaded
regions are obtained by interpolation of the data (the grid is composed of
units of size 0.1 in r.l.u. and 0.02 in meV).

##  VI Theory

In magnetic materials with transition temperatures in the sub-Kelvin range the
dipole-dipole interactions play a significant role. Here we give a brief
theoretical account of the corresponding effects for ytterbium gallium garnet,
for a general discussion see [31, 32]. The Kramers ions Yb3+ in the garnet
structure possess effective $S=1/2$ moments with weakly anisotropic $g$
tensors. Its principal values in the local coordinate frame were obtained in
the EPR experiments on dilute (Y1-xYbx)3Ga5O12 crystals as $g_{x}=2.85$,
$g_{y}=3.60$, and $g_{z}=3.74$ [20]. For simplicity, we adopt an angular
averaged

$g=\sqrt{(g_{x}^{2}+g_{y}^{2}+g_{z}^{2})/3}\approx 3.42\ ,$ (1)

which agrees with the low-temperature magnetization data in Sec. IIC and in
[19]. The garnet structure can be conveniently represented as a bcc lattice
with the primitive unit cell containing 12 magnetic ions residing at

$\bm{\rho}_{1}=\Bigl{(}0,\frac{a}{4},\frac{a}{8}\Bigr{)},\ \
\bm{\rho}_{2}=\Bigl{(}0,\frac{a}{4},\frac{5a}{8}\Bigr{)},\ \
\bm{\rho}_{3,4}=-\bm{\rho}_{1,2}$ (2)

and other eight sites obtained from the above four by cyclic permutations of
their coordinates. The primitive bcc translations are defined as ${\bf
a}_{1}=(-a/2,a/2,a/2)$, ${\bf a}_{2}=(a/2,-a/2,a/2)$, ${\bf
a}_{3}=(a/2,a/2,-a/2)$ with the unit cell volume $V_{\rm uc}=a^{3}/2$. The
lattice constant value $a=12.19$ Å (Sec. V) is used in the following.

###  VI.1 Dipolar interactions and the Curie-Weiss law

The starting point is a standard relation for the magnetic susceptibility
tensor

$\chi^{\alpha\beta}=\frac{(g\mu_{B})^{2}}{k_{B}T}\sum_{i,j}\langle
S^{\alpha}_{i}S^{\beta}_{j}\rangle\ ,$ (3)

expressed via zero-field averages of the angular momentum operators (effective
spins), $k_{B}$ is the Boltzmann constant. We need to consider only the
diagonal elements $\chi^{\alpha\alpha}$, which determine the longitudinal
magnetization.

At high temperatures the magnetic susceptibility can be expanded in powers of
the inverse temperature

$\chi=\frac{C_{1}}{T}+\frac{C_{2}}{T^{2}}+\ldots\simeq\frac{C_{1}}{T-\theta_{CW}}$
(4)

with the Curie constant $C_{1}=(g\mu_{B})^{2}S(S+1)N/3k_{B}$ and the Curie-
Weiss temperature $\theta_{CW}=C_{2}/C_{1}$. The second-order coefficient
$C_{2}$ in the high-temperature series expansion is expressed as

$C_{2}=-\frac{(g\mu_{B})^{2}}{k_{B}^{2}}\sum_{i,j}\Bigl{[}\langle
S^{\alpha}_{i}S^{\alpha}_{j}\hat{\cal H}\rangle_{0}-\langle
S^{\alpha}_{i}S^{\alpha}_{j}\rangle_{0}\langle\hat{\cal H}\rangle_{0}\Bigr{]}\
,$ (5)

where $\hat{\cal H}$ is the magnetic Hamiltonian in the absence of an external
field and $\langle\ldots\rangle_{0}$ denotes paramagnetic averaging. As Eq.
(5) implies, different magnetic interactions contribute additively into
$\theta_{CW}$. Thus, it is possible to separate the dipole shift from other
effects.

The dipole-dipole Hamiltonian is given by

$\hat{\cal H}_{\rm dip}=\frac{1}{2}(g\mu_{B})^{2}\sum_{i,j}\biggl{[}\frac{{\bf
S}_{i}\cdot{\bf S}_{j}}{|{\bf r}_{ij}|^{3}}-\frac{3({\bf S}_{i}\cdot{\bf
r}_{ij})({\bf S}_{j}\cdot{\bf r}_{ij})}{|{\bf r}_{ij}|^{5}}\biggr{]}.$ (6)

Substituting $\hat{\cal H}_{\rm dip}$ into Eq. (5) we obtain the Curie-Weiss
temperature for the component $\chi^{\alpha\alpha}$ as

$\theta^{\alpha\alpha}_{CW}=-\frac{(g\mu_{B})^{2}}{3Nk_{\rm
B}}S(S+1)\sum_{i,j}\biggl{[}\frac{1}{|{\bf
r}_{ij}|^{3}}-\frac{3(r^{\alpha}_{ij})^{2}}{|{\bf r}_{ij}|^{5}}\biggr{]}.$ (7)

Here $N$ is the number of magnetic ions and the sum over both $i$ and $j$ is
kept since the result of a partial summation $\sum_{j}$ may vary for magnetic
ions inside the unit cell. The above expression agrees with the previous works
[33, 34]. We further separate dimensional and dimensionless parts in Eq. (7)
as

$\displaystyle\theta^{\rm dip}_{CW}$ $\displaystyle=$
$\displaystyle-\frac{J_{d}}{3k_{B}}\,S(S+1)\,\frac{a^{3}}{N}\sum_{i,j}\frac{|{\bf
r}_{ij}|^{2}-3z_{ij}^{2}}{|{\bf r}_{ij}|^{5}},$
$\displaystyle\frac{J_{d}}{k_{B}}=\frac{(g\mu_{B})^{2}}{a^{3}}\approx
4.02~{}\textrm{mK}\ .$

The strength of dipolar coupling between adjacent ions is estimated as
$E_{nn}=(8/\sqrt{6})^{3}J_{d}S^{2}/k_{B}=35$ mK.

Lattice sums in Eqs. (7) and (VI.1) are conditionally convergent, i.e., the
outcome depends on the shape of a chosen region. The standard way to treat
this problem is to use the Lorentz sphere construction, which explicitly
separates the shape-dependent demagnetization part [32]. We omit the
demagnetization contribution assuming that a corresponding correction is
included in the experimental data. The rest is represented as a discrete-
lattice sum inside the Lorentz sphere and its surface integral. In cubic
crystals the lattice sum vanishes [this applies only to the double sum in Eq.
(VI.1)] and one is left with the surface integral:

$\frac{1}{N}\sum_{i,j}\frac{|{\bf r}_{ij}|^{2}-3z_{ij}^{2}}{|{\bf
r}_{ij}|^{5}}=-\frac{4\pi}{3v_{0}}\ ,$ (9)

where $v_{0}=a^{3}/24$ is the volume per one magnetic ion in a garnet lattice.
Thus, we obtain for Yb3Ga5O12

$\theta_{CW}^{\rm dip}=32\pi J_{d}\,\frac{S(S+1)}{3k_{B}}\approx 101\
\textrm{mK}\ .$ (10)

This theoretical estimate is very close to the experimental value
$\theta_{CW}^{\rm exp}=97$ mK, Sec. II.3. The remaining small difference can
be attributed to an uncertainty of our experimental measurements, to the Van
Vleck contribution, or an imprecise knowledge of the $g$-factor. Hence, we
conclude that magnetic properties of Yb3Ga5O12 can be well approximated by the
dipolar spin-1/2 model with nearly isotropic $g\sim 3.42$. A possible exchange
coupling between the nearest-neighbor ytterbium ions does not exceed 1–5 % of
the dipole-dipole interaction between them.

###  VI.2 Magnons in the saturated state

In a strong field, magnetic moments become almost parallel to the field
direction. Since we choose an isotropic $g$-tensor for Yb3+ ions, the only
source of a remaining small canting is the dipolar interaction. We neglect
this effect, which is justified once the Zeeman gap $\Delta_{H}$ in the magnon
spectrum is sufficiently large in comparison with $J_{d}$ ($E_{nn}$).

Assuming that the external field $H$ and the local moments are oriented along
the cubic $\hat{\bf z}$ axis we use the Holstein-Primakoff transformation to
transform the Hamiltonian (6) to the quadratic bosonic form. The position of a
magnetic ion is expressed as ${\bf r}_{ni}=\bm{\rho}_{n}+{\bf r}_{i}$, where
$\bm{\rho}_{n}$ is its location in the primitive unit cell and ${\bf r}_{i}$
is the cell position. Then,

$\displaystyle\hat{\cal H}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,J_{d}S\sum_{ni,mj}\Bigl{[}-D^{zz}_{ni,mj}(a_{ni}^{\dagger}a_{ni}+a_{mj}^{\dagger}a_{mj})$
(11) $\displaystyle+$
$\displaystyle\frac{1}{2}(D^{xx}_{ni,mj}+D^{yy}_{ni,mj})(a_{ni}^{\dagger}a_{mj}+a_{mj}^{\dagger}a_{ni})\Bigr{]}$
$\displaystyle+$ $\displaystyle H\sum_{i}a_{ni}^{\dagger}a_{ni}\ .$

where the dipolar matrix is defined as

$D^{\alpha\beta}_{ni,mj}=\frac{\delta^{\alpha\beta}}{|{\bf
r}_{ij}^{nm}|^{3}}-\frac{3(r_{ij}^{nm})^{\alpha}(r_{ij}^{nm})^{\beta}}{|{\bf
r}_{ij}^{nm}|^{5}}$ (12)

with ${\bf r}_{ij}^{nm}={\bf r}_{ni}-{\bf r}_{mj}$ and all distances measured
in units of the lattice constant $a$. We have excluded from Eq. (11) the
anomalous terms $a_{ni}a_{mj}$ ($a^{\dagger}_{ni}a^{\dagger}_{mj})$. These as
well as the linear terms contribute quadratically to the magnon energy
$\varepsilon_{\bf k}$: $\sim J_{d}^{2}/\Delta_{H}$ and can be neglected in
comparison to the linear in $J_{d}$ dispersion originating from the normal
terms.

After the Fourier transformation, the quadratic form $\hat{\cal H}_{2}$
changes to

$\displaystyle\hat{\cal H}_{2}$ $\displaystyle=$ $\displaystyle\sum_{n,{\bf
k}}a_{n{\bf k}}^{\dagger}a_{n\bf
k}\bigl{[}H-J_{d}S\sum_{m}D^{zz}_{nm}(0)\bigr{]}$ $\displaystyle+$
$\displaystyle\frac{1}{2}\,J_{d}S\sum_{n,m,{\bf k}}a_{n{\bf
k}}^{\dagger}a_{m\bf k}\bigl{[}D^{xx}_{nm}({\bf k})+D^{yy}_{nm}({\bf
k})\bigr{]}\ ,$

where the Fourier transform of the dipolar matrix is

$D_{nm}^{\alpha\beta}(\mathbf{k})=\sum_{j}D^{\alpha\beta}_{ni,mj}\;e^{-i\mathbf{k}\cdot({\bf
r}_{j}-{\bf r}_{i})}\ .$ (14)

The matrix elements $D_{nm}^{\alpha\beta}(\mathbf{k})$ are computed using the
Ewald summation technique. The details of this method are described, for
example, in [31]. Below, we give only the final expression in a general form
suitable for any non-Bravais lattice:

$\displaystyle D_{nm}^{\alpha\beta}({\bf k})$ $\displaystyle=$
$\displaystyle-\frac{4k^{3}_{c}}{3\sqrt{\pi}}\,\delta_{\alpha\beta}\delta_{nm}+\frac{4\pi}{V_{\rm
uc}}\sum_{\mathbf{G}}{}^{\prime}\;\frac{(k^{\alpha}+G^{\alpha})(k^{\beta}+G^{\beta})}{|\mathbf{k}+\mathbf{G}|^{2}}\>e^{-|\mathbf{k}+\mathbf{G}|^{2}/4k_{c}^{2}}\>e^{i(\mathbf{k}+\mathbf{G})\cdot(\bm{\rho}_{m}-\bm{\rho}_{n})}$
$\displaystyle+$ $\displaystyle\sum_{\bf R}e^{-i\mathbf{k}\cdot{\bf
R}}\bigg{[}\frac{\textrm{erfc}(k_{c}R_{nm})}{R_{nm}^{3}}\Big{(}\delta_{\alpha\beta}-\frac{3R_{nm}^{\alpha}R_{nm}^{\beta}}{R_{nm}^{2}}\Big{)}-\frac{2k_{c}}{\sqrt{\pi}R_{nm}^{2}}\,e^{-k^{2}_{c}R_{nm}^{2}}\Big{(}2k_{c}^{2}R_{nm}^{\alpha}R_{nm}^{\beta}-\delta_{\alpha\beta}+\frac{3R_{nm}^{\alpha}R_{nm}^{\beta}}{R_{nm}^{2}}\Big{)}\bigg{]}.$

Here, $\textrm{erfc}(x)$ is the complementary error function, the first sum is
performed over the reciprocal lattice points ${\bf G}=m_{1}{\bf
b}_{1}+m_{2}{\bf b}_{2}+m_{3}{\bf b}_{3}$, whereas the second sum is taken
over the direct lattice ${\bf R}=n_{1}{\bf a}_{1}+n_{2}{\bf a}_{2}+n_{3}{\bf
a}_{3}$, ${\bf R}_{nm}={\bf R}+\bm{\rho}_{n}-\bm{\rho}_{m}$, and $k_{c}$ is an
arbitrary cutoff parameter. Independence of a final result on the actual value
of $k_{c}\in[0.2,2]$ can serve as a convenient consistency check. For
arbitrary ${\bf k}\neq{\bf G}$, both sums in Eq. (VI.2) rapidly converge.
However, the term with ${\bf G}=0$ in the first sum becomes ill-determined for
${\bf k}\to 0$ and has to be dropped out, which is indicated by a prime near
the summation symbol. This regularization procedure is equivalent to
neglecting the demagnetization effects, i.e., assuming that a magnetic sample
has cylindrical or slab geometry with a parallelly applied field.

Figure 10: Magnon spectra along the symmetry lines in the Brillouin zone
computed for the dipolar spin model with ${\bf H}\parallel[001]$, $\mu_{0}H=2$
T. The wavevector components are measured in units of $2\pi/a$.

The magnon spectra were computed by numerical diagonalization of Eq. (VI.2).
Figure 10 displays results for an applied field of $\mu_{0}H=2$ T with ${\bf
H}\parallel[001]$ by showing the magnon energy $\varepsilon_{\bf k}$ for the
symmetry directions in the Brillouin zone. Since the primitive unit cell
contains 12 magnetic sites, the excitation spectra have twelve branches with
weak dispersion proportional to $J_{d}$. The branches are gathered into
several groups depending on the momentum direction. The total magnon bandwidth
$\sim 0.04$ meV is an order of magnitude smaller than the gap $\Delta_{H}\sim
0.4$ meV. Such a large separation of the energy scales justifies the
approximation adopted above. However, for weak fields $\mu_{0}H\lesssim 0.4$ T
the neglected anomalous terms as well as the local spin canting have to be
taken into account again.

An interesting feature of the computed spectrum is multiple crossing points
between different magnon branches. Some of them correspond to crossing of
three or even more branches, e.g., at $(1/2,1/2,1/2)$ in Fig. 10(c). Such
conical crossing points may be stable in the cubic garnet crystals, if the
Bloch states in the meeting branches transform according to a three-
dimensional irreducible representation. As a result, magnon bands may possess
interesting topological properties. These questions deserve a separate
theoretical investigation. Comparing with experiment, we note that the average
position of magnon bands at $\sim 0.4$ meV matches well the neutron data of
Sec. V. The limited instrumental resolution as well as thermal fluctuations
($T\sim 2E_{nn}$) cause a smearing effect producing a broad double-peak
structure in the experimental data of Fig. 9 instead of individual magnon
lines. This structure is consistent with the observed tendency of magnon
branches to further bunch into several groups. Still, the experimental
bandwidth is twice larger than the calculated one. Additional future
measurements and calculations may be needed to clarify the origin of the
observed discrepancy.

##  VII Discussion

Previously, the diffuse scattering was measured in the (1 -1 0)–(1 1 2) plane
and analyzed using reverse Monte Carlo calculations performed for the powder
average response [24]. In our work, the diffuse scattering pattern shown in
Fig. 2 is seen projected from a 4-fold axis. Three peaks are found at 0.36,
1.1 and 1.95 Å-1. The first and last values coincide approximately with the
ones found in Ref. 24, which reports another peak at 1.86 Å-1. These typical
ranges of Q are also characteristic of the short range correlations found in
Gd3Ga5O12. This highlights the common spin configurations involved, with
contributions from trimers and 10-spin loops [35]. Furthermore, we show that
the corresponding correlations in Yb3Ga5O12 are significantly increased below
200 mK while a $\bf{Q}$ independent signal (except for magnetic form factor
dependence) is still present at 5 K. The temperature dependence of the elastic
signal finds echo in the increase of the relaxation rate found by $\mu$SR in
the range $T_{\lambda}$ to 0.4 K [23].

The three modes of the magnetic excitation spectrum of Yb3Ga5O12 we observe
are in agreement with results of Sandberg et al. [24] who report rather
similar energies (0.06, 0.12 and 0.7 meV). Our field dependence analysis could
however show that the first mode has a lower energy $E_{0}=0.032$ meV. A
precise determination at zero field would need experiments with higher
resolution. It is to be noted that the elastic and inelastic responses seem
not to be correlated since the strong $\bf{Q}$-dependence of the former is not
reflected in the latter. Interestingly, this magnetic excitation spectrum is
very similar to the one observed in Gd3Ga5O12 where the sequence in energies
is 0.04, 0.12 and 0.58 meV [36], while the two systems are expected to be
different since the Yb3Ga5O12 ground state is supposed to be ordered. This
confirms the similar nature of the correlated state (before ordering) in both
systems.

The intermediate spin configuration we observe under field through the maximum
intensity of the $(1,1,0)$ Bragg peak at 0.15 T presumably corresponding to a
spin canting also reminds the complicated magnetic structure observed in
Gd3Ga5O12 [37, 35]. The studies in this latter compound were nevertheless
hardly conclusive due to the absence of large isotopic single crystals
suitable for neutron scattering experiments. We therefore expect that further
studies in Yb3Ga5O12 will help elucidate the complex field induced magnetic
structure. In contrast to Gd3Ga5O12, we do not observe any signatures of field
induced transitions in magnetization measurements so that the anomalous field
variation of the Bragg peaks could be indicative of a crossover regime. In
addition, our ac susceptibility data do not show evidence for a magnetic
freezing in Yb3Ga5O12 down to 70 mK contrary to the Gd3Ga5O12 case [38].

The most original contribution of the present study is the focus on the spin
dynamics under magnetic field. The experimental ferromagnetic Curie-Weiss
temperature is quantitatively well accounted for the dipolar interactions. The
magnetic excitation spectrum in the saturated phase is also well reproduced
with a magnon calculation based on a dipole-dipole Hamiltonian. A challenge is
the understanding of the spin excitation at zero field and the global magnetic
behavior under field below the saturation. An interesting question is to know
if the dipolar interactions can explain the canting away from the applied
field below the saturated phase as proposed for Gd3Ga5O12 (See Fig. 1(c) of
Ref. 39).

Finally, we have demonstrated that the field dependence of the specific heat
is governed by the field evolution of the band of low-lying magnon
excitations. This highlights presence of an additional microscopic
contribution to the enhanced magnetocaloric effect in Yb3Ga5O12 determined by
collective excitation modes, in line with the theoretical proposal [3].
Indeed, in paramagnet salts that are widely used for ADR applications the
excitation spectrum consists solely of a quasielastic component with width
proportional to $k_{B}T$. Within resolution of the neutron experiments it
appears as an elastic signal [40] and the spectral weight is simply
transferred to a ferromagnetic Bragg peak once the applied field exceeds
$k_{B}T$.

##  VIII Conclusion

Our study demonstrates that the magnetic response of Yb3Ga5O12 is composed of
two distinct contributions: from the ${\bf Q}$-dependent elastic signal and
from almost dispersionless inelastic modes. This landscape appears to be quite
similar to the one found in another garnet, namely Gd3Ga5O12, extensively
studied in the past. The present work reveals the dominant role of dipolar
interactions in Yb3Ga5O12. The theoretical calculations for the spin-1/2
dipolar model reproduce the gross features of the experimental magnon spectrum
in the saturated phase and quantitatively account for the Curie-Weiss
temperature. Thus, ytterbium gallium garnet provides an interesting example of
a quantum dipolar magnet on the hyperkagome lattice. Further experimental and
theoretical investigations of the spin dynamics are necessary to elucidate the
interplay between its highly symmetric frustrated geometry and the long-range
dipolar interactions. In addition, the obtained results confirm that strong
spin correlations in a frustrated geometry can enhance the magnetocaloric
response. We hope that this observation will motivate a continuing search of
suitable frustrated magnetic materials for ADR applications.

##  IX Acknowledgements

The INS data collected at the ILL for the present work are available at Ref.
41, 42. The sample orientation for the neutron experiments was made on IN3 and
Orient Express. E.L., J.-P.B., C.M. and M.E.Z. acknowledge financial support
from ANR, France, Grant No. ANR-18-CE05-0023. J.-P.B., C.M. and M.E.Z.
acknowledge financial support from the Materials transverse program of CEA
“CeramOx”.

## References

  * [1] L. Balents, Nature 464, 199 (2010).
  * [2] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Science 367, eaay0668 (2020).
  * [3] M. E. Zhitomirsky, Phys. Rev. B 67, 104421 (2003).
  * [4] P. Wikus, E. Canavan, S. Trowbridge Heine, K. Matsumoto, and T. Numazawa, Cryogenics 62, 150 (2014).
  * [5] S. S. Sosin, L. A. Prozorova, A. I. Smirnov, A. I. Golov, I. B. Berkutov, O. A. Petrenko, G. Balakrishnan, and M. E. Zhitomirsky, Phys. Rev. B 71, 094413 (2005).
  * [6] M. Orendác, J. Hanko, E. Cizmár, A. Orendácová, M. Shirai, and S. T. Bramwell, Phys. Rev. B 75, 104425 (2007).
  * [7] B. Wolf, U. Tutsch, S. Dörschug, C. Krellner, F. Ritter, W. Assmus, and M. Lang, J. Appl. Phys. 120, 142112 (2016).
  * [8] S. Pakhira, C. Mazumdar, R. Ranganathan, and M. Avdeev, Sci. Reports 7, 7367 (2017).
  * [9] D. Jang, T. Gruner, A. Steppke, K. Mitsumoto, C. Geibel, and M. Brando, Nature Commun. 6, 8680 (2015).
  * [10] Y. Tokiwa, B. Piening, H. S. Jeevan, S. L. Bud’ko, P. C. Canfield, and P. Gegenwart, Sci. Advances 2, e1600835 (2016).
  * [11] D. A. Paixao Brasiliano, J.-M. Duval, C. Marin, E. Bichaud, J.-P. Brison, M. Zhitomirsky, and N. Luchier, Cryogenics 105, 103002 (2020).
  * [12] Y. Tokiwa, S. Bachus, K. Kavita, A. Jesche, A. A. Tsirlin, and P. Gegenwart, Commun Mater 2, 42 (2021).
  * [13] Y. Cai, M. N. Wilson, J. Beare, C. Lygouras, G. Thomas, D. R. Yahne, K. Ross, K. M. Taddei, G. Sala, H. A. Dabkowska, A. A. Aczel, and G. M. Luke, Phys. Rev. B 100, 184415 (2019).
  * [14] R. Wawrzyńczak, B. Tomasello, P. Manuel, D. Khalyavin, M. D. Le, T. Guidi, A. Cervellino, T. Ziman, M. Boehm, G. J. Nilsen, and T. Fennell, Phys. Rev. B 100, 094442 (2019).
  * [15] S. Petit, F. Damay, Q. Berrod, and J. M. Zanotti, Phys. Rev. Research 3, 013030 (2021).
  * [16] I. A. Kibalin, F. Damay, X. Fabréges, A. Gukassov, and S. Petit, Phys. Rev. Research 2, 033509 (2020).
  * [17] J. A. M. Paddison, H. Jacobsen, O. A. Petrenko, M. T. Fernández-Diaz, P. P. Deen, A. L. Goodwin, Science 350, 179 (2015).
  * [18] R. A. Buchanan, K. A. Wickersheim, J. J. Pearson, and G. F. Herrmann, Phys. Rev. 159, 245 (1967).
  * [19] J. Filippi, J. C. Lasjaunias, B. Hebral, J. Rossat-Mignod, and F Tcheou, J. Phys. C: Solid State 13, 1277 (1980).
  * [20] J. W. Carson and R. L. White, J. Appl. Phys. 31, S53 (1960).
  * [21] J. J. Pearson, G. F. Herrmann, K. A. Wickersheim, and R. A. Buchanan, Phys. Rev. 159, 251 (1967).
  * [22] J. A. Hodges, P. Bonville, M. Rams, and K. Królas, J. Phys. Condens. Matter 15, 4631 (2003).
  * [23] P. Dalmas de Réotier, A. Yaouanc, P. C. M. Gubbens, C. T. Kaiser, C. Baines, and P. J. C. King, Phys. Rev. Lett. 91, 167201 (2003).
  * [24] L. O. Sandberg, R. Edberg, K. S. Pedersen, M. C. Hatnean, G. Balakrishnan, L. Mangin-Thro, A. Wildes, B.Fåk, G. Ehlers, G. Sala, P. Henelius, K. Lefmann, P. P. Deen, arXiv:2005.10605.
  * [25] H. W. Capel, Physica 31, 1152 (1965).
  * [26] C. Paulsen, in Introduction to physical techniques in molecular magnetism: Structural and macroscopic techniques (Servicio de Publicaciones de la Universidad de Zaragoza, 2001), p. 1.
  * [27] A. Aharoni, J. Appl. Phys. 83, 3432 (1998).
  * [28] J. R. Stewart, P. P. Deen, K. H. Andersen, H. Schober, J.-F. Barthélémy, J. M. Hillier, A. P. Murani, T. Hayes, and B. Lindenau, J. Appl. Cryst. 42, 69 (2009).
  * [29] The center energy at ${\bf Q}=(\frac{3}{2},\frac{3}{2},0)$ is on average 0.007 meV lower than at ${\bf Q}=(3,0,0)$. For the linear fit, the energy was averaged over the two $\bf{Q}$ positions, and this does not change the results within the accuracy of the measurements.
  * [30] For the field applied along the $c$-axis, (0 1 0) is equivalent to (1 0 0). The former direction was measured due to experimental constraints.
  * [31] K. De’Bell, A. B. MacIsaac, and J. P. Whitehead, Rev. Mod. Phys. 72, 225 (2000).
  * [32] N. Majlis, The Quantum theory of magnetism, 2nd edition (World Scientific Publishing, Singapore, 2007).
  * [33] J. M. Daniels, Proc. Phys. Soc. A 66, 673 (1953).
  * [34] M. J. P. Gingras, B. C. den Hertog, M. Faucher, J. S. Gardner, S. R. Dunsiger, L. J. Chang, B. D. Gaulin, N. P. Raju, and J. E. Greedan, Phys. Rev. B 62, 6496 (2000).
  * [35] P. P. Deen, O. Florea, E. Lhotel, and H. Jacobsen, Phys. Rev. B 91, 014419 (2015).
  * [36] P. P. Deen, O. A. Petrenko, G. Balakrishnan, B. D. Rainford, C. Ritter, L. Capogna, H. Mutka, and T. Fennell, Phys. Rev. B 82, 174408 (2010).
  * [37] O. A. Petrenko, G. Balakrishnan, D. McK Paul, M. Yethiraj, G. J. McIntyre, and A. S. Wills, J. Phys.: Conf. Series 145, 012026 (2009).
  * [38] P. Schiffer, A. P. Ramirez, D. A. Huse, P. L. Gammel, U. Yaron, D. J. Bishop and A. J. Valentino, Phys. Rev. Lett. 74, 2379 (1995).
  * [39] N. d’Ambrumenil, O. A. Petrenko, H. Mutka, and P. P. Deen, Phys. Rev. Lett. 114, 227203 (2015).
  * [40] A. T. Boothroyd, Principles of Neutron Scattering from Condensed Matter, p. 194 (Oxford University Press, Oxford, 2020).
  * [41] https://doi.ill.fr/10.5291/ILL-DATA.5-42-501
  * [42] https://doi.ill.fr/10.5291/ILL-DATA.4-05-724

##  Appendix A Fit of the Neutron data

The elastic line does not evolve between 0.25 and 2 T and therefore all the
magnetic elastic scattering is suppressed at 0.25 T. The remaining incoherent
scattering is then estimated from the 2 T data and is described by two
gaussians of width 0.039 meV and 0.096 meV, the second one having an intensity
of 1/18 of the first one (This reflects globally the non-gaussian profiles of
the monochromator and analyser transfer functions). In all the fits of the
elastic line, only an overall intensity (with the second gaussian intensity
fixed at 1/18 of the first one) and an offset are free parameters. The
energies of the modes given in the text are corrected for the offset. For the
inelastic modes under magnetic field, the intensity and peak positions are
always free parameters. The width of the lowest energy mode is fixed at 0 T to
the value found at 0.25 T, where the peak is resolved and the width is then a
free parameter for all finite fields. For the second mode, the width is fixed
at 0, 0.25 and 0.5 T and free at 1 T when it merges with the third mode. For
the third mode, all parameters are free.


# A constrained cosmological model in $f(R,L_{m})$ gravity

J. K. Singh<EMAIL_ADDRESS>Department of Mathematics, Netaji Subhas
University of Technology, New Delhi-110078, India  Shaily
<EMAIL_ADDRESS>Department of Mathematics, Netaji Subhas University of
Technology, New Delhi-110078, India  Ratbay Myrzakulov<EMAIL_ADDRESS>Ratbay Myrzakulov Eurasian International Centre for Theoretical Physics,
Astana, Kazakhstan.  Harshna Balhara<EMAIL_ADDRESS>Department of
Mathematics, Netaji Subhas University of Technology, New Delhi-110078, India

###### Abstract

In this article, we study the expanding nature of universe in the contest of
$f(R,L_{m})$ gravity theory, here $R$ represents the Ricci scalar and $L_{m}$
is the matter Lagrangian density. With a specific form of $f(R,L_{m})$, we
obtain the field equations for flat FLRW metric. We parametrize the
deceleration parameter in terms of the Hubble parameter and from here we find
four free parameters, which are constraints and estimated by using $H(z)$,
$Pantheon$, and their joint data sets. Further, we investigate the evolution
of the deceleration parameter which depicts a transition from the deceleration
to acceleration phases of the universe. The evolution behaviour of energy
density, pressure, and EoS parameters shows that the present model is an
accelerated quintessence dark energy model. To compare our model with the
$\Lambda$CDM model we use some of the diagnostic techniques. Thus, we find
that our model in $f(R,L_{m})$ gravity supports the recent standard
observational studies and delineates the late-time cosmic acceleration.

Keywords: $f(R,L_{m})$ gravity, Observational constraints, FLRW metric, dark
energy

##  I Introduction

The main objective of current and future cosmological investigations is to
better understand the nature of cosmic acceleration by putting the cold dark
matter ($\Lambda$CDM) standard cosmological model and its variations to the
test [1, 2, 3, 4]. Despite the fact that the cosmological constant $\Lambda$
matches the data from observations well, it has two main flaws: the
coincidence problem and the cosmological constant problem [5]. One may create
gravity theories that are equal to $\Lambda$CDM at the background level but
display unique and exciting signs on perturbation dynamics by using the
cosmological constant as the underlying source of the accelerated behaviour.
In light of this idea, we will investigate if there is a gravity theory that
possesses these qualities and can prove to be better than the $\Lambda$CDM
paradigm.

One of the most straightforward options is to include an auxiliary function of
the Ricci scalar $R$ in the action, which results in the $f(R)$ theory [6, 7,
8]. $f(R)$ theories of gravity are distinct among higher-order gravity
theories in that they appear to be the only ones that can prevent the well-
known and disastrous Ostrogradski instability [9, 10, 11]. When one quantizes
the fluctuations of the scalar field in the background metric, one obtains the
heavy scalar particles as well as the graviton in the case of $f(R)$ gravity.
Given the size of the scalar particles in $f(R)$ gravity, the pressure may be
negligible and the strength of the interaction between these scalar particles
and the ordinary matter should be on par with that of gravitational forces.
Such a scalar particle may therefore be a prime candidate for dark matter.
Furthermore, because $f(R)$ theories do not contain ghosts, they can be
selected so that the additional degrees of freedom compared to those of GR do
not inescapably lead to serious viability issues. Additionally, it has been
seen that some models display chameleon behaviour, enabling the theory to have
cosmic effects that explain for the universe’s current acceleration [13, 12,
14, 15]. Severe weak field restrictions in the Solar System regime appear to
rule out the majority of the models put forth so far [16, 17, 18, 19, 20],
while workable models do exist[21, 22, 23]. Along with the solar system and
equivalence principle restrictions on $f(R)$ gravity, observational signals of
$f(R)$ dark energy models have been provided in [24, 25, 26, 27]. It has been
discussed in [28, 29, 30] that there are more $f(R)$ models that unify early
inflation and dark energy and pass local tests. To learn more about the many
cosmological consequences of $f(R)$ gravity models, see references [31, 32,
33].

The explicit coupling of an arbitrary function of the Ricci scalar R with the
matter Lagrangian density $L_{m}$ was included in the theory as a
generalization of the $f(R)$ gravity theories [34]. This model was expanded to
include the scenario of arbitrary couplings in both matter and geometry [35].
The non-minimal matter-geometry coupling’s effects on cosmology and
astrophysics have been thoroughly researched in [36, 37, 38, 39, 40]. Of all
the gravitational theories built in Riemann space, the $f(R,L_{m})$ gravity
theory can be viewed as their maximum extension [41]. In this gravity theory,
test particle motion is non-geodesic and an additional force orthogonal to the
four velocity vectors appears. Gonclaves and Moraes considered the
$f(R,L_{m})$ gravity when analysing cosmology from non-minimal matter geometry
coupling [42]. Wu et al. studied Constraints of energy conditions and DK
instability criterion on $f(R,L_{m})$ gravity models [43]. Wang and Liao have
recently investigated the energy conditions in $f(R,L_{m})$ gravity [44].

The structure of the article is as follows: In Section II, we derive Einstein
Fields equations for $f(R,L_{m})$ theory of gravity and take the parametrized
deceleration parameter in terms of the Hubble parameter. In the next section,
we constrain the best fit values of free parameters by using Hubble data set
$H(z)$ (77 points), Pantheon data set (1048 points) and joint data set
$(H(z)+Pantheon)$. In Section IV, we investigate the evolution profile of
various cosmological parameters and diagnostics techniques. Finally, in the
last section, we summarize the outcomes of our work.

##  II Einstein Field Equations

The gravitational interactions in $f(R,L_{m})$ modified theory of gravity can
be written as [41]

$A=\int f(R,L_{m})\sqrt{-g}d^{4}x$ (1)

where $f(R,L_{m})$ is a random function of the Ricci scalar $R$ and the matter
Lagrangian density term $L_{m}$. By taking the variation in the $L_{m}$ term
we can obtain the energy momentum tensor of the matter and the outcome is

$T_{ij}=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}L_{m})}{\delta
g^{ij}}=g_{ij}L_{m}-2\frac{\partial L_{m}}{\partial g^{ij}},$ (2)

By taking the variation in action, Eq. (1) leads the expression as

$\delta A=\int\Bigg{[}f^{R}(R,L_{m})\delta R+f^{L_{m}}(R,L_{m})\frac{\delta
L_{m}}{\delta g^{ij}}\delta g^{ij}-\frac{1}{2}g_{ij}f(R,L_{m})\delta
g^{ij}\Bigg{]}\sqrt{-g}d^{4}x,$ (3)

where $f^{R}=\frac{\partial f}{\partial R}$ and $f^{L_{m}}=\frac{\partial
f}{\partial L_{m}}$. Variation in Ricci scalar can be calculated from the
following manner

$\delta R=\delta(g^{ij}R_{ij})=R_{ij}\delta
g^{ij}+g^{ij}(\nabla_{k}\delta\Gamma^{k}_{ij}-\nabla_{j}\delta\Gamma^{k}_{ik}),$
(4)

where $\nabla_{k}$ is the covariant w.r.t $g_{ij}$ and using the variated
value of Christoffel symbol i.e.
$\delta\Gamma^{k}_{ij}=\frac{1}{2}g^{k\zeta}(\nabla_{i}\delta
g_{j\zeta}+\nabla_{j}\delta g_{\zeta i}-\nabla_{\zeta}\delta g_{ij})$, Eq. (4)
leads to

$\delta R=R_{ij}\delta g^{ij}+g_{ij}\nabla_{i}\nabla^{i}\delta
g^{ij}-\nabla_{i}\nabla_{j}\delta g^{ij}.$ (5)

Therefore, from Eq. (3) we obtain the Einstein field equations for
$f(R,L_{m})$ gravity as

$f^{R}(R,L_{m})R_{ij}+(g_{ij}\nabla_{i}\nabla^{i}-\nabla_{i}\nabla_{j})f^{R}(R,L_{m})-\frac{1}{2}[f(R,L_{m})-f^{L_{m}}(R,L_{m})L_{m}]g_{ij}=\frac{1}{2}f^{L_{m}}(R,L_{m})T_{ij}.$
(6)

Here, the energy momentum tensor for perfect fluid filled universe is written
as $T_{ij}=(\rho+p)u_{i}u_{j}+pg_{ij}$, where p is used for isotropic
pressure, $\rho$ is matter energy density and the component of four velocities
i.e. $u_{i}=(1,0,0,0)$. We use the flat FLRW metric line element to study the
present model and it can be expressed as

$ds^{2}=-dt^{2}+a^{2}(t)(dx^{2}+dy^{2}+dz^{2}),$ (7)

we obtain the EFEs from Eq. (6) as follows

$3H^{2}f^{R}+\frac{1}{2}(f-f^{R}R-f^{L_{m}}L_{m})+3H\dot{f}^{R}=\frac{1}{2}f^{L_{m}}\rho,$
(8)

and

$\dot{H}f^{R}+3H^{2}f^{R}-\ddot{f}^{R}-3H\dot{f}^{R}+\frac{1}{2}(f^{L_{m}}L_{m}-f)=\frac{1}{2}f^{L_{m}}p.$
(9)

In $f(R,L_{m})$ theory of gravity many models have been discussed and the
motivation behind the investigation of our model is the related work of Harko
et al. [45]. Taking this work into account, we consider the following
functional form as:

$f(R,L_{m})=\frac{R}{2}+L_{m}^{\lambda}+\zeta,$ (10)

where $\lambda$ and $\zeta$ are arbitrary constants.

The Einstein field Eqs. (8) and (9) in $f(R,L_{m})$ gravity takes the form as:

$3H^{2}=\rho^{\lambda}(2\lambda-1)-\zeta$ (11)
$4\dot{H}+9H^{2}=n\rho^{n-1}p+\beta-\rho^{n}(n-1)$ (12)

As the above two field equations have three unknowns namely $H$, $\rho$, and
$p$, therefore we need an additional constraint to get the solution. In this
regard we have a crucial dimensionless element in comprehending the expansion
scenario of the universe, that is deceleration parameter $q$. The negative
value of $q$ indicates the accelerating universe and the positive value of $q$
shows the deceleration phase of the universe. If $q=0$ then it shows a
constant rate of expansion of the universe. Hence, with this motivation, we
express the deceleration parameter in terms of the Hubble parameter $H$ as
[46]

$q=\alpha-\frac{\beta}{H}$ (13)

where $\alpha$ and $\beta$ are model parameters with $\beta>0$. The scale
factor $a$ is solved by using (13) and the relation
$q=-\frac{a\ddot{a}}{\dot{a}^{2}}$ as

$a(t)=n(e^{\beta t}-1)^{\frac{1}{1+\alpha}}$ (14)

where $n$ is the free parameter. Also Eq. (14) leads to the following solution
for the Hubble parameter

$H=\frac{\beta e^{\beta t}}{(1+\alpha)(e^{\beta t}-1)}$ (15)

Since the scale factor $a$ in terms of redshift $z$ is expressed as
$\frac{a}{a_{0}}=\frac{1}{1+z}$ where $a_{0}=1$ is the present value. Thus,
using this relation Hubble parameter in terms of redshift can be written as

$H=\frac{\beta\left((n(z+1))^{\alpha+1}+1\right)}{\alpha+1}.$ (16)

And so Eqs. (11), (12) and (16) leads the values of $\rho$ and $p$ as

$\rho=\left(\frac{(\alpha+1)^{2}\zeta+3\left(\beta+\beta(n(z+1))^{\alpha+1}\right)^{2}}{(\alpha+1)^{2}(2\lambda-1)}\right)^{1/\lambda},$
(17)

and

##  III Statistical observation for model parameters

In this section, we present the observational data and modality to calculate
the free model parameters. To find the best fit value of model parameters, we
use Hubble data set $H(z)$ (77 points), Pantheon data set (1048 points) and
joint data set $(H(z)+Pantheon)$. For this methodology, we can obtain the
value of Hubble parameter from Eq. (16) in terms of present age $(t_{0})$ and
$H_{0}$ as

$H=H_{0}(1-e^{-\beta t_{0}})[{n(1+z)}^{1+\alpha}+1],$ (18)

Where $H_{0}=67.4km/s/Mpc$ [4] and we constraint remaining model parameters
$n$, $\alpha$, $\beta$ and $t_{0}$ using Markov Chain Monte Carlo (MCMC)
method with emcee library in python for the mentioned data sets.

###  III.1 Hubble Data set

In the present model, we have taken 77 points data set from $H(z)$ assessment.
$H(z)$ data is obtained from cosmic chronometers and very useful to examine
the dark section of the universe [47]. By minimizing the value of chi-square,
we can get the best fit results for the model parameters. Therefore,
$\chi^{2}$ is written as

$\chi_{Hub}^{2}(n,\alpha,\beta,t_{0})=\sum\limits_{i=1}^{77}\frac{[H(n,\alpha,\beta,t_{0},z_{i})-H_{obs}(z_{i})]^{2}}{\sigma_{z_{i}}^{2}},$
(19)

where $H(n,\alpha,\beta,t_{0},z_{i})$ and $H_{obs}$ stands for the theoretical
and observed values of Hubble parameter respectively and $\sigma_{(z_{i}}$
indicates the standard deviation for each $H(z_{i})$.

(a)  Figure 1: Likelihood contours for $H(z)$ data set. Table 1: Constrained values of model parameters. Dataset |  $n$ |  $\alpha$ |  $\beta$ |  $t_{0}$
---|---|---|---|---
$H(z)$ |  $0.4165^{+0.0036}_{-0.0035}$ |  $1.637^{+0.056}_{-0.056}$ |  $1.51^{+0.29}_{-0.29}$ |  $13.00^{+0.55}_{-0.55}$
$Pantheon$ |  $2.1802^{+0.0013}_{-0.0013}$ |  $2.57^{+0.85}_{-0.85}$ |  $3.9994^{+0.0006}_{-0.0002}$ |  $13.12^{+0.82}_{-0.54}$
H(z)+Pantheon |  $0.4999^{+0.0001}_{-0.0001}$ |  $1.3019^{+0.0001}_{-0.0001}$ |  $2.7284^{+0.0001}_{-0.0001}$ |  $13.1006^{+0.0001}_{-0.0001}$

###  III.2 Pantheon Data

In furtherance of study observational data and to constrain the free
parameters namely $n$, $\alpha$, $\beta$ and $t_{0}$, we use the currently
accessible $SNeIa$ sample, which is known as Pantheon sample [48]. In total
Pantheon data contains 1048 points, which is the outcome of compilation of
various $SNeIa$ surveys. The CfA1-CfA4 surveys, the Carnegie Supernova Project
(CSP), the Sloan Digital Sky Survey (SDSS), the Supernovae Legacy Survey
(SNLS), the Pan-STARRS1 (PS1) and many Hubble Space Telescope (HST) samples
incorporate in it [49, 50, 51, 52, 53]. In all these surveys, the range of
redshift lies between 0.01 to 2.26. Supernovae data mainly use to investigate
the expansion rate of the universe, thus for this probe we introduce the
apparent magnitude as

$m(z)=M+5log_{10}\Bigg{[}\frac{d_{L}(z)}{1Mpc}\Bigg{]}+25,$ (20)

where $M$ is the absolute magnitude and $D_{L}$ is the luminosity distance.
For flat universe $D_{L}$ can be obtain by

$d_{L}(z)=(1+z)\int_{0}^{z}\frac{c}{H(z^{*})}dz^{*},$ (21)

where $c$ is the speed of light. Also, the apparent magnitude can be expressed
in terms of Hubble free luminosity distance
$(D_{L}(z)=\frac{H_{0}d_{L}(z)}{c})$ as

$m(z)=M+5log_{10}[D_{L}(z)]+5log_{10}\Bigg{[}\frac{c/H_{0})}{Mpc}\Bigg{]}+25.$
(22)

Here a degeneracy between $M$ and $H_{0}$ can be notice and since distance
modulus $\mu=m-M$, therefore $\chi^{2}$ for Pantheon data can be written in
terms of $\mu$ as

$\chi_{SN}^{2}(n,\alpha,\beta,t_{0})=\sum\limits_{i=1}^{1048}\left[\frac{\mu_{th}(n,\alpha,\beta,t_{0},z_{i})-\mu_{obs}(z_{i})}{\sigma_{\mu(z_{i})}}\right]^{2},$
(23)

where $\mu_{th}$ and $\mu_{obs}$ used for theoretical and observed distance
modulus respectively. Here, we perform Markov chain Monte Carlo method to
constrain the best fit value of free parameters, which are written in Table-I.

(a)  Figure 2: Likelihood contours for $Pantheon$ data set.

###  III.3 Joint Data sets ($H(z)+Pantheon$)

For joint analysis, the expression for $\chi^{2}$ function is given as

$\chi_{HS}^{2}(n,\alpha,\beta,t_{0})=\chi_{Hub}^{2}(n,\alpha,\beta,t_{0})+\chi_{SN}^{2}(n,\alpha,\beta,t_{0}).$
(24)

(a)  Figure 3: Likelihood contours for $H(z)+Pantheon$ data set.

Joint data analysis is useful to get the stronger constraint value of the free
parameters and the obtained values of $n$, $\alpha$, $\beta$ and $t_{0}$, from
all the observational data sets are placed in Table-I.

##  IV Dynamical nature of the Universe

In this section, we scrutinize various cosmological parameters with the
constraint values of free parameters obtained from the Hubble data set,
Pantheon data set and their combined data set. To understand the evolution of
the universe in a better way, we study cosmic parameters one by one. From Eqs.
(13) and (16), we can calculate the value of the deceleration parameter as

$q=\alpha-\frac{\alpha+1}{(n(z+1))^{\alpha+1}+1}.$ (25)

From Fig. 4(a), we observe that the deceleration parameter changes its phase
from positive to negative i.e. universe transit from deceleration to
acceleration. It is noticed that for the Hubble data set and joint data set,
the universe was decelerating in the past and accelerating at present and in
the future too. But the Pantheon sample shows deceleration in the past and
present and afterwards, it accelerates in late times, which does not support
recent evidence which suggests that the Universe’s expansion is speeding up.

(a)  (b)

Figure 4: The graphs of $q$ and $j$ vs. $z$.

Now, to explore the nature of the universe in the present model, we study the
parameters which contain the higher-order derivatives of $a$. Here we discuss
the jerk parameter ($j$), which is also known as a jolt, pulse, bounce,
impulse, surge etc. [54]. The value of $j$ can be calculated from the formula
as

$j=\frac{\dddot{a}}{aH^{3}},$ (26)

and using Eqs. (14) and (16), jerk parameter can be obtain in terms of
redshift as

$j=\frac{\alpha\left(\alpha+(2\alpha+1)(n(z+1))^{\alpha+1}-1\right)(n(z+1))^{\alpha+1}+1}{\left((n(z+1))^{\alpha+1}+1\right)^{2}}.$
(27)

Fig. 4(b) indicates that from early times to late times the value of $j$
decreases for all observation samples and finally approaches to 1, which shows
that this model is different from $\Lambda CDM$ model at the early Universe
and similar to $\Lambda CDM$ model in late times.

To understand the physical behaviour of the universe we study various
parameters like energy density $(\rho)$, isotropic pressure $(p)$, equation of
state parameter $(\omega)$ etc. Using the expressions given in Eqs. (17) and
(LABEL:18), the plots of energy density and pressure are drawn in Fig. 5 and
here we fixed the value of $\lambda$ and $\zeta$ as $0.55$ and $20$
respectively. Where the first figure shows that for all taken data sets,
initially energy density is very high and afterwards it is receding from early
times to late times. Fig. 5(b) depicts that for this model pressure is
negative throughout the range of redshift for the mentioned data sets and
according to the standard cosmology, negative pressure indicates the
accelerated expansion of the universe.

(a)  (b)  (c)

Figure 5: The graphs of $\rho_{\phi}$, $p_{\phi}$ and $\omega$ vs. $z$ for
distinct observational data.

The value of the EoS parameter can be obtained from the equation
$\omega=\frac{p}{\rho}$ in terms of $z$, where the value of $p$ and $\rho$ are
taken from Eqs. (17) and (LABEL:18). In Fig. 5(c), trajectories of
$\omega~{}vs.~{}z$ depict that for $H(z)$ and joint data sets, the value of
$\omega$ lies in the quintessence region for the whole range of $z$ and for
$Pantheon$ data set our model is in the phantom region at early times as well
as at present but in future, it enters in quintessence region and remains
there.

Now another tool to compare our model with $\Lambda$CDM model and to
understand the geometric behavior of the universe for the present model is
statefinder diagnostic techniques [55]. In this technique, two geometrical
diagnostic parameters $(r,s)$ are introduced and formulated as:

$r=\frac{\dddot{a}}{aH^{3}},$ (28)

$s=\frac{r-1}{3(q-\frac{1}{2})},~{}~{}~{}where~{}~{}q\neq\frac{1}{2}$ (29)

(a)  (b)

Figure 6: The plots of $s-r$ and $q-r$.

Statefinder technique is useful to compare various dark energy model with
$\Lambda$CDM, $SCDM$, $HDE$ model. Now with the help of $s-r$ trajectories it
is quite easy to understand the nature of our model. In $s-r$ plot, the point
$(0,1)$ represents the $\Lambda$CDM model, $(2/3,1)$ denotes holographic dark
energy $(HDE)$ model and $(1,1)$ represents standard cold dark matter $(SCDM)$
model. In our model, we plot the $s-r$ curve by calculating the values of $s$
and $r$ from Eqs. (14), (16) and (25). We notice that all the $s-r$
trajectories are converging to $\Lambda$CDM. As the arrows are indicating that
trajectories for each observational data enters from Chaplygin Gas region and
finally approach to $\Lambda$CDM model. The present values of $(s,r)$ are
$(0.0079,0.9703)$, $(0.0424,0.8589)$, $(2.3920,14.6384)$ for $H(z)$,
$Pantheon$ and their joint data sets respectively, which shows that $H(z)$,
$Pantheon$ and joint data trajectories are very close to $\Lambda$CDM at
present (see Fig. 6(a)).

In Fig. 6(b), we examine the $q-r$ curves for each datasets, where the point
$(1,\frac{1}{2})$ denotes the $SCDM$ model, $(-1,1)$ represents steady state
$(SS)$ model and the dotted horizontal line $r=1$ indicates $\Lambda$CDM
model. In our model, we observe that all the trajectories converge to $SS$
model after crossing the transition line ($\Lambda$CDM line). Thus, our model
is similar to steady state model in late times as the evolution of the
trajectories approaches to $SS$ in future.

The Om diagnostic technique is generally used to observe the variations of
various dark energy models from $\Lambda$CDM model and [56, 57]. This tool can
be discussed without using the EoS parameter and is defined in terms of $z$
and $H(z)$ as:

$Om=\frac{\Bigg{(}\frac{H(z)}{H_{0}}\Bigg{)}^{2}}{Z^{3}+3z^{2}+3z}.$ (30)

According to modern cosmology, positive curvature of the trajectory indicates
the Phantom model, negative curvature of the trajectory shows the quintessence
model and if there is no curvature then it indicates $\Lambda$CDM model. In
the present model, Fig. 7 shows that our model is a quintessence dark energy
model for a large range of $z$.

(a)  Figure 7: The plot of Om diagnostic  vs. $z$.

##  V Conclusion

In this article, we have studied the late-time behavior of the universe for a
flat isotropic FLRW metric in $f(R,L_{m})$ theory of gravity. Now, to find a
suitable solution set of the Einstein Field equations, we parametrize the
deceleration parameter in terms of Hubble parameter $H$ and evaluate the scale
factor $a$ as a function of $t$. We constrained the free parameters using
Hubble $H(z)$, $Pantheon$, and their joint data sets. The best-fitted values
of the free parameters are placed in Table-I. In addition, we have
investigated the various other cosmological parameters using the observational
values. At present, deceleration parameter shows the acceleration phase of the
Universe for $H(z)$ and joint data sets but for $Pantheon$ data picture is
little different (see Fig. 4(a)). At late times our model behaves like
$\Lambda$CDM and this is the outcome of jerk parameter as the value of j
tending to $1$ for all taken data sets.

In this model, the energy density decreases monotonically as $z$ decreases and
at late times it diminishes to zero at late time. Isotropic pressure increased
from high redshift to low redshift, and it is negative throughout the range of
$z$, which indicates the accelerating expansion of the universe. Also, the
evolution profile of the EoS parameter $\omega$ shows that the curves of
$\omega$ lie in the quintessence region for $H(z)$ and joint data sets, and it
$\omega$ changes its value from phantom to quintessence region for $Pantheon$
data sets but all the three curves show the quintessence dark energy model at
late time. Furthermore, we investigated statefinder diagnostics where the
$s-r$ plot also shows that at late times our model converges to $\Lambda$CDM
model and $q-r$ trajectory represents that at the end our model converges to
$SS$ model. Finally, the evolution profile of the Om diagnostic indicates that
our is a quintessence model.

## References

  * [1] A. G. Riess et al. [Supernova Search Team], Astron. J. 116, 1009-1038 (1998).
  * [2] S. Perlmutter et al. [Supernova Cosmology Project], Astrophys. J. 517, 565-586 (1999).
  * [3] P. A. R. Ade et al. [Planck], Astron. Astrophys. 594, A13 (2016).
  * [4] N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)].
  * [5] E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753-1936 (2006).
  * [6] H. A. Buchdahl, Mon. Not. Roy. Astron. Soc. 150, 1 (1970).
  * [7] H. Kleinert and H. J. Schmidt, Gen. Rel. Grav. 34, 1295-1318 (2002).
  * [8] R. Kerner, Gen. Rel. Grav. 14, 453-469 (1982).
  * [9] T. P. Sotiriou, J. Phys. Conf. Ser. 189, 012039 (2009).
  * [10] A. de la Cruz-Dombriz and A. Dobado, Phys. Rev. D 74, 087501 (2006).
  * [11] T. H. Loo, A. De, S. Arora and P. K. Sahoo, Eur. Phys. J. C 82, 705 (2022).
  * [12] S. Capozziello, S. Nojiri, S. D. Odintsov and A. Troisi, Phys. Lett. B 639, 135-143 (2006).
  * [13] S. M. Carroll, V. Duvvuri, M. Trodden and M. S. Turner, Phys. Rev. D 70, 043528 (2004).
  * [14] L. Amendola, D. Polarski and S. Tsujikawa, Phys. Rev. Lett. 98, 131302 (2007).
  * [15] S. Nojiri and S. D. Odintsov, Phys. Rev. D 74, 086005 (2006).
  * [16] T. Chiba, Phys. Lett. B 575, 1-3 (2003).
  * [17] A. L. Erickcek, T. L. Smith and M. Kamionkowski, Phys. Rev. D 74, 121501 (2006).
  * [18] T. Chiba, T. L. Smith and A. L. Erickcek, Phys. Rev. D 75, 124014 (2007).
  * [19] S. Nojiri and S. D. Odintsov, Phys. Lett. B 659, 821-826 (2008).
  * [20] G. J. Olmo, Phys. Rev. D 75, 023511 (2007).
  * [21] W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007).
  * [22] I. Sawicki and W. Hu, Phys. Rev. D 75, 127502 (2007).
  * [23] L. Amendola and S. Tsujikawa, Phys. Lett. B 660, 125-132 (2008).
  * [24] S. Tsujikawa, Phys. Rev. D 77, 023507 (2008).
  * [25] T. Liu, X. Zhang and W. Zhao, Phys. Lett. B 777, 286-293 (2018).
  * [26] S. Capozziello and S. Tsujikawa, Phys. Rev. D 77, 107501 (2008).
  * [27] A. A. Starobinsky, JETP Lett. 86, 157-163 (2007).
  * [28] S. Nojiri and S. D. Odintsov, Phys. Lett. B 657, 238-245 (2007).
  * [29] S. Nojiri and S. D. Odintsov, Phys. Rev. D 77, 026007 (2008).
  * [30] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani and S. Zerbini, Phys. Rev. D 77, 046009 (2008).
  * [31] J. Santos, J. S. Alcaniz, M. J. Reboucas and F. C. Carvalho, Phys. Rev. D 76, 083513 (2007).
  * [32] S. Capozziello, V. F. Cardone and V. Salzano, Phys. Rev. D 78, 063504 (2008).
  * [33] R. C. Nunes, S. Pan, E. N. Saridakis and E. M. C. Abreu, JCAP 01 (2017), 005.
  * [34] O. Bertolami, C. G. Boehmer, T. Harko and F. S. N. Lobo, Phys. Rev. D 75, 104016 (2007).
  * [35] T. Harko, Phys. Lett. B 669, 376-379 (2008).
  * [36] T. Harko, Phys. Rev. D 81, 084050 (2010).
  * [37] T. Harko, Phys. Rev. D 81, 044021 (2010).
  * [38] S. Nesseris, Phys. Rev. D 79, 044015 (2009).
  * [39] V. Faraoni, Phys. Rev. D 76, 127501 (2007).
  * [40] V. Faraoni, Phys. Rev. D 80, 124040 (2009).
  * [41] T. Harko and F. S. N. Lobo, Eur. Phys. J. C 70, 373-379 (2010).
  * [42] B. S. Gonçalves and P. H. R. S. Moraes, [arXiv:2101.05918 [gr-qc]].
  * [43] Y. B. Wu, Y. Y. Zhao, Y. Y. Jin, L. L. Lin, J. B. Lu and X. Zhang, Mod. Phys. Lett. A 29, 1450089 (2014).
  * [44] J. Wang and K. Liao, Class. Quant. Grav. 29, 215016 (2012).
  * [45] T. Harko and F. S. N. Lobo, Galaxies 2, no. 3, 410-465 (2014).
  * [46] R. K. Tiwari, D. Sofuoğlu and V. K. Dubey, Int. J. Geom. Meth. Mod. Phys. 17 (2020) no.12, 2050187.
  * [47] L. Chimento and M. I. Forte, Phys. Lett. B 666, 205-211 (2008).
  * [48] D. M. Scolnic et al. [Pan-STARRS1], Astrophys. J. 859, no.2, 101 (2018).
  * [49] A. G. Riess, R. P. Kirshner, B. P. Schmidt, S. Jha, P. Challis, P. M. Garnavich, A. A. Esin, C. Carpenter, R. Grashius and R. E. Schild, et al. Astron. J. 117, 707-724 (1999).
  * [50] S. Jha, R. P. Kirshner, P. Challis, P. M. Garnavich, T. Matheson, A. M. Soderberg, G. J. M. Graves, M. Hicken, J. F. Alves and H. G. Arce, et al. Astron. J. 131, 527-554 (2006).
  * [51] M. Hicken, P. Challis, S. Jha, R. P. Kirsher, T. Matheson, M. Modjaz, A. Rest and W. M. Wood-Vasey, Astrophys. J. 700, 331-357 (2009).
  * [52] C. Contreras, M. Hamuy, M. M. Phillips, G. Folatelli, N. B. Suntzeff, S. E. Persson, M. Stritzinger, L. Boldt, S. Gonzalez and W. Krzeminski, et al. Astron. J. 139, 519-539 (2010).
  * [53] M. Sako et al. [SDSS], Publ. Astron. Soc. Pac. 130, no.988, 064002 (2018).
  * [54] J. K. Singh, A. Singh, G. K. Goswami and J. Jena, Annals Phys. 443, 168958 (2022).
  * [55] U. Alam, V. Sahni, T. D. Saini and A. A. Starobinsky, Mon. Not. Roy. Astron. Soc. 344, 1057 (2003).
  * [56] V. Sahni, A. Shafieloo and A. A. Starobinsky, Phys. Rev. D 78, 103502 (2008).
  * [57] C. Zunckel and C. Clarkson, Phys. Rev. Lett. 101, 181301 (2008).