# Running Hubble constant
from the SNe Ia Pantheon sample?

Tiziano Schiavone  Department of Physics “E. Fermi”, University of Pisa, Polo
Fibonacci, Largo B. Pontecorvo 3, I-56127 Pisa, Italy
<EMAIL_ADDRESS>INFN, Istituto Nazionale di Fisica Nucleare,
Sezione di Pisa, Polo Fibonacci, Largo B. Pontecorvo 3, I-56127, Pisa, Italy
Giovanni Montani  ENEA, Fusion and Nuclear Safety Department, C.R. Frascati,
Via E. Fermi 45, I-00044 Frascati (Roma), Italy  Physics Department,
“Sapienza” University of Rome, P.le Aldo Moro 5, I-00185 Roma, Italy  Maria
Giovanna Dainotti  National Astronomical Observatory of Japan, 2 Chome-21-1
Osawa, Mitaka, Tokyo 181-8588, Japan
<EMAIL_ADDRESS>The Graduate University for Advanced Studies,
SOKENDAI, Shonankokusaimura, Hayama, Miura District, Kanagawa 240-0193, Japan
Space Science Institute, Boulder, CO, USA  Biagio De Simone  Department of
Physics “E.R. Caianiello”, University of Salerno, Via Giovanni Paolo II, 132,
Fisciano, I-84084 Salerno, Italy  INFN Gruppo Collegato di Salerno, Sezione di
Napoli, c/o Dipartimento di Fisica “E.R. Caianiello”, Ed. F, Università di
Salerno, Via Giovanni Paolo II, 132, Fisciano, I-84084 Salerno, Italy  Enrico
Rinaldi  Physics Department, University of Michigan, Ann Arbor, MI 48109, USA
Theoretical Quantum Physics Laboratory, Center for Pioneering Research, RIKEN,
2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Interdisciplinary Theoretical &
Mathematical Science Program, RIKEN (iTHEMS), 2-1 Hirosawa, Wako, Saitama,
351-0198, Japan  Gaetano Lambiase  Department of Physics “E.R. Caianiello”,
University of Salerno, Via Giovanni Paolo II, I-132-84084 - Fisciano, Salerno
- Italy  INFN Gruppo Collegato di Salerno, Sezione di Napoli, c/o Dipartimento
di Fisica “E.R. Caianiello”, Ed. F, Università di Salerno, Via Giovanni Paolo
II, 132, Fisciano, I-84084 Salerno, Italy

###### Abstract

The mismatch between different independent measurements of the expansion rate
of the Universe is known as the Hubble constant ($H_{0}$) tension, and it is a
serious and pressing problem in cosmology. We investigate this tension
considering the dataset from the Pantheon sample, a collection of 1048 Type Ia
Supernovae (SNe Ia) with a redshift range $0<z<2.26$. We perform a binned
analysis in redshift to study if the $H_{0}$ tension also occurs in SNe Ia
data. Hence, we build equally populated subsamples in three and four bins, and
we estimate $H_{0}$ in each bin considering the $\Lambda$CDM and
$w_{0}w_{a}$CDM cosmological models. We perform a statistical analysis via a
Markov Chain Monte Carlo (MCMC) method for each bin. We observe that $H_{0}$
evolves with the redshift, using a fit function
$H_{0}(z)=\tilde{H}_{0}(1+z)^{-\alpha}$ with two fitting parameters $\alpha$
and $\tilde{H}_{0}$. Our results show a decreasing behavior of $H_{0}$ with
$\alpha\sim 10^{-2}$ and a consistency with no evolution between 1.2 $\sigma$
and 2.0 $\sigma$. Considering the $H_{0}$ tension, we extrapolate $H_{0}(z)$
until the redshift of the last scattering surface, $z=1100$, obtaining values
of $H_{0}$ consistent in 1 $\sigma$ with the cosmic microwave background (CMB)
measurements by Planck. Finally, we discuss possible $f(R)$ modified gravity
models to explain a running Hubble constant with the redshift, and we infer
the form of the scalar field potential in the dynamically equivalent Jordan
frame.

##  I Introduction

The Hubble constant tension is one of the biggest open problems of modern
cosmology. Several independent measurements of the actual expansion rate of
the Universe, the Hubble constant $H_{0}$, provide incompatible values. The
status of all these discrepancies and theoretical attempts is summarized in a
review [1]. It should be noted that this inconsistency occurs between
measurements referred to the early Universe, and those based on local probes
(for instance, Cepheids or Type Ia supernovae - SNe Ia) in the late Universe.
This tension has become more and more serious, little by little more precise
measurements have reduced gradually the error bars. More precisely, there is a
discrepancy in 4.4 $\sigma$ between the Planck data [2], $H_{0}=67.4\pm
0.5\,\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}$, and the value obtained from
Cepheids [3], $H_{0}=74.03\pm 1.42\,\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}$.
This tension could point out inconsistencies in Planck data or local probes
or, alternatively, it could be a signal for a new cosmology.

Regarding the latter point, modified gravity theories are studied to try to
solve open problems in cosmology, such as the Hubble constant tension or the
nature of dark energy, which cannot be fully explained in General Relativity
(GR) and the standard $\Lambda$CDM cosmological model. In particular, the
$f(R)$ modified gravity theories [4, 5] provide a generalization of GR,
including an extra geometrical degree of freedom given by a function $f$ of
the Ricci scalar $R$, which implies a cosmological dynamics that differs from
GR. The scalar-tensor theories in the so-called Jordan frame [4, 5, 6] are
dynamically equivalent to the $f(R)$ proposals; within this framework, the
additional scalar degree of freedom is provided by a scalar field, which is
non-minimally coupled to the metric.

Since the $H_{0}$ tension concerns measurements referred to the early and late
Universe, in this work we investigate the possibility of a hidden evolutionary
effect that could imply a running Hubble constant with the redshift [7, 8, 9,
10], focusing on the data analysis of the Pantheon sample [11]. In this
regard, we follow a binned analysis of the Pantheon sample in three and four
redshift bins to extract the values of $H_{0}$ in each bin. Then, we observe
an unexpected decreasing trend of $H_{0}$ with the redshift, and we interpret
this extra degree of freedom via a modified gravity scenario.

This work is organized as follows: in Sec. 1 we recall briefly some basic
notions about cosmology using SNe Ia datasets; in Sec. 2 we introduce the
$f(R)$ extended models in the Jordan frame; in Sec. 3 we follow a binning
approach to address the Hubble constant tension within the redshift range of
the Pantheon sample; in Sec. 4 we discuss our results and the theoretical
interpretations in the $f(R)$ modified cosmology.

We adopt the metric signature $\left(-,+,+,+\right)$, and we set the speed of
light $c=1$. The Einstein constant is denoted with $\chi\equiv 8\,\pi\,G$, and
$G$ is the Newton constant.

##  II Type Ia Supernova cosmology

SNe Ia are usually regarded as standard candles since they are characterized
by their uniformity in the absolute magnitude profiles. Hence, they are very
useful to estimate cosmic distances. In this regard, we recall the definition
of the theoretical distance modulus of a SN:

$\mu_{\textrm{th}}=m_{th}-M=5\hskip 2.15277ptlog_{10}\ d_{L}(z)+25\,,$ (1)

where $m_{th}$ and $M$ are the apparent and absolute magnitude of a SN,
respectively, and $d_{L}$ is the luminosity distance, which depends on a
specific cosmological model. In a flat geometry, $d_{L}$ is given by [12]

$d_{L}\left(z\right)=\left(1+z\right)\,\int_{0}^{z}\frac{dz^{\prime}}{H\left(z^{\prime}\right)}\,.$
(2)

Considering a dark energy component in which the equation of state parameter
$w$ evolves with the redshift $z$, i.e. $w=w(z)$, and neglecting relativistic
components in the late Universe, the Hubble function $H(z)$ contained in Eq.
(2) is written as

$H\left(z\right)=H_{0}\,\sqrt{\Omega_{m0}\,\left(1+z\right)^{3}+\Omega_{DE\,0}\,\exp\left[3\,\int_{0}^{z}\left[1+w\left(z^{\prime}\right)\right]\frac{dz^{\prime}}{1+z^{\prime}}\right]}\,,$
(3)

where $\Omega_{m0}$, $\Omega_{DE\,0}\equiv 1-\Omega_{m0}$ are the cosmological
density parameters of the matter and dark energy components, respectively, at
the present redshift $z=0$. Note that if $w=-1$, a cosmological constant
$\Lambda$ is reproduced, and the standard $\Lambda$CDM scenario is recovered.
Differently, considering in Eq. (3) the $w_{0}w_{a}$CDM model, given by the
Chevallier-Polarski-Linder (CPL) parameterization [13, 14]
$w\left(z\right)=w_{0}+w_{a}\times z/\left(1+z\right)$ with parameters $w_{0}$
and $w_{a}$, we end up in a slowly evolving dark energy scenario slightly
different from a cosmological constant.

Moreover, it should be noted that the Hubble constant $H_{0}$ is degenerate
with $M$, as you can see easily combining Eqs. (1) (2) and (3). Specifically,
$M$ is calibrated to $-19.35$ in the Pantheon sample such that
$H_{0}=70.0\,\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}$ [11].

To discriminate between several cosmological models, it is useful to compare
$\mu_{\text{th}}$ in Eq. (1) with the observed distance modulus
$\mu_{\text{obs}}$, which is obtained from SNe Ia data,

$\mu_{\textrm{obs}}=m_{B}-M+\alpha\,x_{1}-\beta\,c+\Delta M+\Delta B\,,$ (4)

where $m_{B}$ is the B-band apparent magnitude, and $M$ is the absolute
magnitude in the B-band$x_{1}$ of a SN with $c=0$ and $x_{1}=0$. SNe
parameters are the color $c$ and the stretch $x_{1}$, while $\alpha$ and
$\beta$ are coefficients. Finally, $\Delta M$ is a distance correction related
to the host-galaxy mass, and $\Delta B$ is a bias correction. The Pantheon
sample data is available in the repository by Scolnic et al. (2018) [11]
(https://github.com/dscolnic/Pantheon).

To perform a statistical analysis, we define
$\Delta\mu=\mu_{\textrm{obs}}-\mu_{\textrm{th}}$ as the difference between the
distance moduli in Eqs. (1) and (4). Furthermore, the $\chi^{2}$ is built as

$\chi^{2}=\Delta\mu^{T}\,C^{-1}\,\Delta\mu\,,$ (5)

where $C=C_{sys}+D_{stat}$ is the full covariance matrix. The $C_{sys}$ matrix
involves systematic errors, while $D_{stat}$ is a diagonal matrix that
contains statistical errors for each SN, due to peculiar velocities,
photometry, bias, lensing, and intrinsic scatter [11]. Referring to the
Pantheon sample, $\Delta\mu$ represents a vector with $1048$ components, while
$C$, $C_{sys}$ and $D_{stat}$ are $1048\times 1048$ square matrices.

##  III f(R) modified cosmology in the Jordan frame

Within the framework of the $f(R)$ modified gravity [4, 5], the gravitational
Lagrangian density becomes a function $f(R)$ of the scalar curvature $R$,
involving an additional scalar degree of freedom with respect to GR. Adopting
the metric formalism, it can be noted that the generalized field equations are
of fourth-order in the metric [5], and if $f\left(R\right)=R$ one might obtain
the Einstein-Hilbert equations in GR. The scalar-tensor theories in the so-
called Jordan frame [4, 5, 6] are useful to rewrite $f(R)$ models in a
dynamically equivalent form, which provides field equations of lower order.
The total action in the Jordan frame is written as

$S_{J}=\frac{1}{2\,\chi}\,\int
d^{4}x\,\sqrt{-g}\,\left[\phi\,R-V\left(\phi\right)\right]+S_{M}\left(g_{\mu\nu},\psi\right)\,,$
(6)

where $g$ is the determinant of $g_{\mu\nu}$ the metric tensor, $S_{M}$ is the
matter contribution to the total action, and $\psi$ is referred to the matter
fields. Moreover, in this paradigm, the extra degree of freedom given by
$f\left(R\right)$ is expressed as a scalar field $\phi$, defined as
$\phi=f^{\prime}\left(R\right)=df/dR$, which is non-minimally coupled to the
metric, and it is governed by the scalar field potential
$V\left(\phi\right)=\phi\,R\left(\phi\right)-f\left(R\left(\phi\right)\right)$.

The field equations in the Jordan frame for a flat FLRW geometry in the late
Universe, considering a pressureless matter component and neglecting
relativistic species, are given by the following equation system [6]:

$\displaystyle
H^{2}=\frac{\chi\,\rho}{3\,\phi}-H\,\frac{\dot{\phi}}{\phi}+\frac{V\left(\phi\right)}{6\,\phi}$
(7a)
$\displaystyle\frac{\ddot{a}}{a}=-\frac{\chi\,\rho}{6\,\phi}+\frac{V\left(\phi\right)}{6\,\phi}-\frac{H}{2}\,\frac{\dot{\phi}}{\phi}-\frac{1}{2}\,\frac{\ddot{\phi}}{\phi}$
(7b) $\displaystyle
3\ddot{\phi}-2\,V\left(\phi\right)+\phi\,\frac{dV}{d\phi}+9\,H\,\dot{\phi}=\chi\,\rho,$
(7c)

where $\rho\left(t\right)$ is the matter energy density, and a dot denotes a
time derivative. The first equation above is the modified Friedmann equation,
the second one is the modified acceleration equation, while the third equation
describes the scalar field dynamics. Note that all these equations are now of
second order in the equivalent formalism in the Jordan frame, but the
disadvantage is the introduction of a non-minimally coupling between the
scalar field and the metric.

It should be emphasized that the gravitational coupling constant is restated
in the Jordan frame. More specifically, looking at Eq. (7a),
$G_{\text{eff}}=G/\phi$ may be regarded as an effective gravitational
coupling. Modified gravity models are very interesting to explore alternative
scenarios for the cosmic acceleration in the late Universe without a true
cosmological constant (see the models by [5, 15, 16, 17]).

Finally, following the Hu-Sawicki formalism [15], we introduce auxiliary
variables $y_{H}$ and $y_{R}$ that are useful to write the luminosity distance
in the $f(R)$ gravity. We write the Hubble parameter and Ricci scalar as

$H^{2}=m^{2}\left(\left(1+z\right)^{3}+y_{H}\right),\hskip
43.05542ptR=m^{2}\left(3\left(1+z\right)^{3}+y_{R}\right),$ (8)

where the two dimensionless variables $y_{H}$ and $y_{R}$ encompass extra
contributions with respect to the matter component in the $\Lambda$CDM model.
We have defined $m^{2}=\chi\,\rho_{0}/3=H_{0}^{2}/\Omega_{m0}$ with $\rho_{0}$
present matter density. Note that if $y_{H}$ is simply a constant, the first
relation is nothing more the Friedmann equation in the $\Lambda$CDM model. In
the context of $f(R)$ gravity, instead, $y_{H}$ and $y_{R}$ evolve and their
dynamics strongly depend on the form of the $f(R)$ function or, equivalently,
the scalar field potential $V\left(\phi\right)$. It is possible to rewrite the
generalized Friedmann equation (7a) and the Ricci scalar in terms of $y_{H}$
and $y_{R}$ and their derivatives to obtain a fist-order differential equation
system (see [15]), which can be solved numerically. Then, we can write the
luminosity distance in a $f(R)$ model

$d_{L}(z)=\frac{(1+z)}{H_{0}}\int^{z}_{0}\frac{dz^{\prime}}{\sqrt{\Omega_{m0}\biggr{(}y_{H}(z^{\prime})+(1+z^{\prime})^{3}\biggr{)}}}$
(9)

in terms of the deviation $y_{H}(z)$, combining Eqs. (2) and (8).

##  IV Binned analysis

We decide to follow a binned approach using the Pantheon sample dataset to
evaluate if the Hubble constant evolves with the redshift [9, 10], motivated
by the fact the independent measurements of $H_{0}$ referred to different
redshifts had shown a tension.

Therefore, we have split the redshift range of the Pantheon sample in bins
each having the same number $N$ of SNe Ia. More precisely, considering three
redshift bins, we have $N\approx 349$ SNe for each bin, and the redshift
ranges for three bins are $0.01<z<0.18$, $0.18<z<0.34$, and $0.34<z<2.26$. In
the case of four bins, instead, $N=262$ in each bin and the redshift ranges
are: $0.01<z<0.13$, $0.13<z<0.25$, $0.25<z<0.42$, and $0.42<z<2.26$. Then, we
have built the subvectors $\Delta\mu$ with $N$ SNe, according to the binning
division and redshift order. We have also divided the full covariance matrix
$C$ into $N\times N$ submatrices for each redshift bin. In this regard,
concerning simply the statistical contribution $D_{stat}$ in $C$, each
diagonal element of the $D_{stat}$ matrix is related to a single SN, hence we
can quickly build the submatrices starting from $D_{stat}$. However, if we
include the systematic errors, we need to consider also $C_{sys}$, which is
not a diagonal matrix. Thus, we used a customized code to select $C_{sys}$
elements involved only with the SNe within the redshift bin taken into
account. Hence, we built properly the submatrices, considering both
statistical and systematic errors related to the SNe for each bin.

To focus only on $H_{0}$ in a one-dimensional analysis and constrain it in
each bin, we decide to fix the value of the other cosmological parameters.
Thus, we assume the fiducial value [11] $\Omega_{m0}=0.298$ for a flat
$\Lambda$CDM model, and $\Omega_{m0}=0.308$ with $w_{0}=-1.009$ and
$w_{a}=-0.129$ for a flat $w_{0}w_{a}$CDM model for all the redshift ranges.
Moreover, we need to fix the absolute magnitude $M$ of SNe Ia in each bin.
Bearing in mind the degeneracy between $M$ and $H_{0}$, we decide to calibrate
$M$ to obtain a value of $H_{0}=73.5\,\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}$
in the first redshift bins, according to measurements of the local probes at
low redshifts.

Thus, we perform a preliminary analysis for both $\Lambda$CDM and
$w_{0}w_{a}$CDM models to obtain the respective value of $M$ for each binning
division in the first redshift bin, which is characterized by a mean redshift
$z\ll 1$. To do this, we minimize the $\chi^{2}$ in Eq. (5). Then, we use the
MCMC methods to sample a posterior distribution. Finally, $M$ is obtained for
each binning division (see Table 1). After this preliminary analysis, we fix
$M$ in the other bins. One might expect that also $H_{0}$ assumes the same
values in all the redshift bins, due to its degeneracy with $M$.

Flat $\Lambda$CDM Model, Fixed $\Omega_{m0}=0.298$
---
Bins | $\tilde{H}_{0}$ | $\alpha$ | $\frac{\alpha}{\sigma_{\alpha}}$ | $M$ | $H_{0}\left(z=11.09\right)$ | $H_{0}\left(z=1100\right)$
| $\,\left(\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}\right)$ |  |  |  | $\,\left(\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}\right)$ | $\,\left(\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}\right)$
3 | $73.577\pm 0.106$ | $0.009\pm 0.004$ | $2.0$ | $-19.245\pm 0.006$ | $72.000\pm 0.805$ | $69.219\pm 2.159$
4 | $73.493\pm 0.144$ | $0.008\pm 0.006$ | $1.5$ | $-19.246\pm 0.008$ | $71.962\pm 1.049$ | $69.271\pm 2.815$
Flat $w_{0}w_{a}$CDM Model, Fixed $\Omega_{m0}=0.308$, $w_{0}=-1.009$ and
$w_{a}=-0.129$
Bins | $\tilde{H}_{0}$ | $\alpha$ | $\frac{\alpha}{\sigma_{\alpha}}$ | $M$ | $H_{0}\left(z=11.09\right)$ | $H_{0}\left(z=1100\right)$
| $\,\left(\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}\right)$ |  |  |  | $\,\left(\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}\right)$ | $\,\left(\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}\right)$
3 | $73.576\pm 0.105$ | $0.008\pm 0.004$ | $1.9$ | $-19.244\pm 0.005$ | $72.104\pm 0.766$ | $69.516\pm 2.060$
4 | $73.513\pm 0.142$ | $0.008\pm 0.006$ | $1.2$ | $-19.246\pm 0.004$ | $71.975\pm 1.020$ | $69.272\pm 2.737$
Table 1: Fitting parameters ($\tilde{H}_{0}$ and $\alpha$) of $H_{0}(z)$ (Eq.
(10) and extrapolated values at high redshifts, after a binned analysis of the
Pantheon sample, focusing on a flat $\Lambda$CDM model (upper part) and a flat
$w_{0}w_{a}$CDM model (lower part). The compatibility of $\alpha$ with zero
(no evolution) is expressed in terms of 1 $\sigma$ in the fourth column. The
new absolute magnitude $M$, which provides $H_{0}=73.5\,\textrm{km
s}^{-1}\,\textrm{Mpc}^{-1}$ in the first bins, is listed in the fifth column.
Finally, the extrapolated values of $H_{0}(z)$ at the redshift of the most
distant galaxies, $z=11.09$, and the last scattering surface, $z=1100$, are
listed in the last two columns. All the errors are presented in 1 $\sigma$.

Figure 1: Evolution of $H_{0}(z)$ with redshift according to Eq. (10) as a
result of a binned analysis of the Pantheon sample for a flat $\Lambda$CDM
model (two upper plots) and a flat $w_{0}w_{a}$CDM model (two lower plots).

Now, we can approach the main part to check the values of $H_{0}$ after a
binned analysis of the Pantheon Sample in three and four bins for both the
$\Lambda$CDM and $w_{0}w_{a}$CDM models. We use again the minimization of
$\chi^{2}$, as illustrated in the preliminary analysis, and in the MCMC we set
the priors: $60$ $\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}$ $<H_{0}<80$
$\textrm{km s}^{-1}\,\textrm{Mpc}^{-1}$. Finally, we have extracted the values
of $H_{0}$ for the all bins. To investigate an evolution of $H_{0}$ with the
redshift $z$, we apply a non-linear fit [9] of $H_{0}(z)$ written as

$H_{0}(z)=\frac{\tilde{H}_{0}}{(1+z)^{\alpha}}\,,$ (10)

with two fitting parameters, $\tilde{H}_{0}$ and $\alpha$. Note that
$\tilde{H}_{0}=H_{0}\left(z=0\right)$. Moreover, if $\alpha\neq 0$, we have an
evolutionary trend, otherwise $H_{0}$ is a constant.

The fitting parameters and their respective errors in 1 $\sigma$ are listed in
Table 1. Note that, for instance, the $\alpha$ parameter is consistent with
zero, namely no evolution, in 2.0 $\sigma$ in the $\Lambda$CDM model using
three bins. In Fig. 1 we can observe a slowly and unexpected decreasing trend
[9] for $H_{0}\left(z\right)$ in the $\Lambda$CDM model and also in the
$w_{0}w_{a}$CDM model. It should be noted that all the $\alpha$ coefficients
are mutually compatible in 1 $\sigma$, therefore pointing out a reliable
decreasing trend of $H_{0}$ with the redshift in the Pantheon sample.

If the observed trend of $H_{0}(z)$ is intrinsic and does not depend on a
specific sample, one may ask what happens with other local probes or in the
early Universe at very high redshifts. Hence, we extrapolate the fit function
of $H_{0}(z)$ at the redshift of the most distant galaxies, $z=11.09$ [18],
and at the redshift of the last scattering surface, $z=1100$, to compare the
latter values of $H_{0}(z)$ with the one inferred from Planck measurements for
the CMB. We find that the extrapolated values (see Table 1) are consistent
within 1 $\sigma$ with the value of $H_{0}$ obtained with the Planck
measurements for both the $\Lambda$CDM and $w_{0}w_{a}$CDM models regardless
the number of bins.

##  V Interpretation of the results in the f(R) modified cosmology

The evolution of the Hubble constant with the redshift needs a physical
interpretation. In this section, we do not focus on possible astrophysical
reasons (see [9]), but we investigate our results in modified gravity
theories.

The Hubble constant must be a constant by definition, and its observed
evolution could be a signal of a wrong framework, for instance, a possible
hidden function of the redshift, which has not been taken into account so far.
As already discussed in [19, 9], a varying Einstein constant $\chi=8\pi G$ (or
a varying Newton constant $G$) may in principle lead to an evolution of
$H_{0}\left(z\right)$. If we look at the function of $H_{0}\left(z\right)$ in
Eq. (10), we require an effective Einstein constant
$\chi\sim\left(1+z\right)^{-2\,\alpha}$ to preserve constant the present
critical density $\rho_{c0}=3\,H_{0}^{2}/\chi$.

We recall that an effective Einstein constant is naturally obtained in the
Jordan frame (Sec. 2) through the rescaling: $\chi\rightarrow\chi\,/\,\phi$.
The extra degree of freedom in the Jordan frame, i.e. the presence of a non-
minimally coupled scalar field $\phi$, which depends on the redshift, provides
indeed an evolving Einstein constant. Moreover, another reason to focus on the
Jordan frame to account for the observed trend of $H_{0}(z)$ lies in the
modified Friedmann equation (7a). Indeed, the Hubble parameter $H(z)$ is
related to matter density $\rho$ via an effective Einstein constant, hence the
dynamics of the scalar field $\phi$ may imply an evolution of the Hubble
constant with redshift. More specifically, we require the following ansatz

$\phi\left(z\right)=\left(1+z\right)^{2\,\alpha}$ (11)

to account for $H_{0}\left(z\right)$ in Eq. (10). These concepts suggest that
an extra degree of freedom with respect to GR might imply the unexpected trend
of $H_{0}\left(z\right)$.

We do not mention so far which scalar field potential we should have in the
Jordan frame, mimicking a cosmological constant in a slow-roll and at the same
time providing an effective Hubble constant. Considering the cosmological
dynamics in the Jordan frame, we can infer $V\left(\phi\right)$ [10], if we
assume the behavior of $H_{0}\left(z\right)$ in Eq. (10) and the ansatz above
for $\phi\left(z\right)$ in Eq. (11). Therefore, the modified Friedmann
equation (7a) allows us to write $V\left(\phi\right)$ as

$V\left(\phi\right)=6\,(1-2\alpha)\,\left(\frac{dz}{dt}\right)^{2}\,\phi^{1-1/\alpha}-6m^{2}\,\phi^{3/2\alpha}\,,$
(12)

where we have used the standard definition of redshift, and
$\rho\sim(1+z)^{3}$ for a matter component. Now, we need to specify the term
$dz/dt$, which is given by

$\frac{dz}{dt}=-(1+z)\,H(z).$ (13)

Here, we can not simply use the extended Friedmann equation (7a) to replace
$H(z)$, because we should set a specific scalar field potential to solve
$H(z)$. Then, we decide to impose the Hubble parameter

$H\left(z\right)=\frac{\tilde{H}_{0}}{(1+z)^{\alpha}}\sqrt{\Omega_{m0}\,\left(1+z\right)^{3}+1-\Omega_{m0}}\,,$
(14)

as suggested from our binned analysis in Sec. 3, although its analytical form
is different from the Hubble parameter obtained using the modified Friedmann
equation. In other words, to obtain $V\left(\phi\right)$, we impose in the
Jordan frame dynamics the same physical effect observed from the redshift
binned analysis of the Pantheon sample.

Figure 2: Scalar field potential $V(\phi)$ in the Jordan frame given by Eq.
(15) and inferred from the evolution of $H_{0}\left(z\right)$ with redshift
(10). The quantity $V(\phi)/m^{2}$ is dimensionless. Note the presence of a
flat region of the potential for $0<z\lesssim 0.3$ or $\phi\lesssim 1.005$.

Hence, combining Eqs. (12), (13) and (14), we obtain the final expression of
the potential [10]:

$\frac{V\left(\phi\right)}{m^{2}}=6\,(1-2\alpha)\,\frac{1-\Omega_{m0}}{\Omega_{m0}}-12\,\alpha\,\phi^{\frac{3}{2\alpha}}\,,$
(15)

where we have considered the relation $\Omega_{m0}=m^{2}/\tilde{H}_{0}^{2}$.
You can see the form of the scalar field potential in Fig. 2, where we fixed
$\alpha=0.009$ (see Table 1) and $\Omega_{m0}=0.298$ for plotting
$V\left(\phi\right)/m^{2}$. Note the occurrence of a flat region of
$V\left(\phi\right)$ for $0<z\lesssim 0.3$, mimicking a cosmological constant
in the era dominated by dark energy, thus validating our approach and the
inferred scalar field potential.

Furthermore, we can obtain the $f(R)$ function related to $V\left(\phi\right)$
in Eq. (15). Using the field equation $R=dV/d\phi$ in the Jordan frame, and
also the relation $f(R)=R\,\phi(R)-V\left(\phi(R)\right)$ [5], we write:

$f\left(R\right)=-6\,m^{2}\,\left[\left(3-2\alpha\right)\,\left(-\frac{R}{18\,m^{2}}\right)^{\frac{3}{3-2\alpha}}+\left(1-2\alpha\right)\frac{1-\Omega_{0m}}{\Omega_{0m}}\right]\,.$
(16)

Expanding this expression for $\alpha\sim 0$, according to the values of
$\alpha$ from the binned analysis (see Table 1), we obtain [10]:

$f\left(R\right)\approx\left(R-6\,m^{2}\,\frac{1-\Omega_{0m}}{\Omega_{0m}}\right)+\frac{2}{3}\alpha\,\left[R\,\ln{\left(-\frac{R}{m^{2}}\right)}-\left(1+\ln{18}\right)\,R+18m^{2}\,\frac{1-\Omega_{0m}}{\Omega_{0m}}\right]+O\left(\alpha^{2}\right)\,.$
(17)

It should be emphasized that the first term above is precisely the
gravitational Lagrangian density in GR with a cosmological constant
$\Lambda=3m^{2}\left(1-\Omega_{0m}\right)\,/\,\Omega_{0m}$, where we used
$m^{2}=\tilde{H}^{2}_{0}\,\Omega_{0m}$. Moreover, it is clear that the first-
order term in $\alpha$ in Eq. (17) gives the deviation from GR. Hence, the
$\alpha$ coefficient stresses corrections to GR.

The expressions inferred for the scalar field potential in Eq. (15), and the
corresponding $f(R)$ function in Eq. (17), are viable in the late Universe
since in all the calculations we do not include relativistic components.
However, it could be interesting to test this form of $V\left(\phi\right)$
using other probes in the late Universe.

The proposed discussions above show a possible interpretation of our results
coming from the binned analysis, and we suggest a candidate for
$V\left(\phi\right)$, assuming the evolution of $H_{0}(z)$.

##  VI Conclusions

In this work, performing a redshift binned analysis of the Pantheon sample of
SNe Ia, we have seen that there is a slow evolution of the Hubble constant
$H_{0}$ with the redshift $z$, described by Eq. (10) with a parameter $\alpha$
that is consistent with zero (no evolution) between 1.2 $\sigma$ and 2.0
$\sigma$. The evolutionary behavior of $H_{0}(z)$ occurs regardless of the
binning division both in the $\Lambda$CDM and $w_{0}w_{a}$CDM models. It seems
that the Hubble constant tension also emerges locally in the narrow redshift
range $0<z<2.26$ of the Pantheon sample of SNe Ia.

Then, extrapolating the fit function of $H_{0}(z)$ (10) at higher redshifts,
we obtain values of $H_{0}$ that are remarkably consistent within 1 $\sigma$
with the CMB measurements by Planck. Our results could point out an intrinsic
trend of $H_{0}(z)$, which might imply the mismatch between independent
measurements of the Hubble constant referred to the early and late Universe.

We remark that we perform a one-dimensional analysis to obtain constraints
only on $H_{0}$, following [9]. Actually, in [10] we extend the analysis, and
we obtain similar results, considering more than one variable in the MCMC
methods, for instance $H_{0}$ and $\Omega_{m0}$, to avoid fixing all the
cosmological parameters except $H_{0}$.

Our analysis could suggest the presence of a hidden astrophysical effect,
which is the reason for the evolution of $H_{0}(z)$ and has not been taken
into account so far. Alternatively, the observed effect may point out that we
need a cosmology beyond the $\Lambda$CDM paradigm. We have discussed the
$f(R)$ modified gravity theories, and we suggest a possible form for the
scalar field potential (15) in the Jordan frame, and the respective $f(R)$
function (17). Furthermore, we mention that a new binned analysis of the
Pantheon sample could be worthwhile, using a specific $f(R)$ model and the
modified version of the luminosity distance (9) to test new physics.

Concerning the future perspectives, to explore further the possible evolution
of cosmological parameters, it could be very useful to include in our analysis
also probes at high redshifts, such as the Gamma-Ray Bursts (GRBs), which in
principle could allow to enlarge the redshift range of the actual Hubble
diagram [20, 21, 22, 23, 24, 25]. This work could be regarded as a preliminary
study for the application of GRBs in cosmology, using, for instance, the
prompt-afterglow relations with the plateau emission [26, 27, 28, 29, 30, 31,
32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42], similarly to the magnetar model
[43, 44, 45].

In conclusion, we have presented a new perspective to address the Hubble
constant tension through a binning approach, thus we foster investigations on
the astrophysical parameters of SNe Ia and further theoretical discussions.

## References

  * [1] E. Di Valentino, O. Mena, P. Supriya, et al., In the realm of the Hubble tension — a review of solutions, Class. Quantum Grav. 38, 153001 (2021)
  * [2] N. Aghanim, Y. Akrami, M. Ashdown, et al., Planck 2018 results. VI. Cosmological parameters, A$\&$A 641, A6 (2020)
  * [3] M. J. Reid, D. W. Pesce and A. G. Riess, An Improved Distance to NGC 4258 and Its Implications for the Hubble Constant, ApJL 886, L27 (2019)
  * [4] S. Nojiri, and S. D. Odintsov, Introduction to modified gravity and gravitational alternative for dark energy, IJGMM 4, 115 (2007)
  * [5] T. P. Sotiriou, and V. Faraoni, $f\left(R\right)$ theories of gravity, Rev. Mod. Phys. 82, 451 (2010)
  * [6] S. Capozziello, and V. Faraoni, Beyond Einstein gravity (Springer, 2013)
  * [7] C. Krishnan, E. Ó Colgáin, Ruchika, et al., Is there an early Universe solution to Hubble tension?, Phys. Rev. D 102, 103525 (2020)
  * [8] C. Krishnan, E. Ó Colgáin, M.M. Sheikh-Jabbari, and Tao Yang, Running Hubble tension and a H0 diagnostic, Phys. Rev. D 103, 103509 (2021)
  * [9] M. G. Dainotti, B. De Simone, T. Schiavone, et al., On the Hubble Constant Tension in the SNe Ia Pantheon Sample, ApJ 912(2), 150 (2021)
  * [10] M. G. Dainotti, B. De Simone, T. Schiavone, et al., On the Evolution of the Hubble constant with the SNe Ia Pantheon Sample and Baryon Acoustic Oscillations: a Feasibility Study for GRB-Cosmology in 2030, Galaxies 10, 24 (2022)
  * [11] D. M. Scolnic, D. O. Jones, A. Rest, et al., The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample, ApJ 859, 101 (2018)
  * [12] S. Weinberg, Cosmology (Oxford University Press, 2008)
  * [13] M. Chevallier, and D. Polarski, Accelerating Universes with Scaling Dark Matter, IJMPD 10, 213 (2001)
  * [14] E. V. Linder, Exploring the Expansion History of the Universe, Phys. Rev. Lett. 90, 091301 (2003)
  * [15] W. Hu, and I. Sawicki, Models of $f(R)$ cosmic acceleration that evade solar system tests, Phys. Rev. D 76, 064004 (2007)
  * [16] A. A. Starobinsky, Disappearing cosmological constant in $f(R)$ gravity, JETP Lett. 86, 157 (2007)
  * [17] S. Tsujikawa, Observational signatures of $f(R)$ dark energy models that satisfy cosmological and local gravity constraints, Phys. Rev. D 77, 023507 (2008)
  * [18] P. A. Oesch, G. Brammer, P. G. van Dokkum, et al., A Remarkably Luminous Galaxy at z=11.1 Measured with Hubble Space Telescope Grism Spectroscopy, ApJ 819, 129 (2016)
  * [19] L. Kazantzidis, and L. Perivolaropoulos, Hints of a local matter underdensity or modified gravity in the low $z$ Pantheon data, Phys. Rev. D 102, 023520 (2020)
  * [20] V. F. Cardone, S. Capozziello, and M. G. Dainotti, An updated gamma-ray bursts Hubble diagram, MNRAS 400, 775 (2009)
  * [21] V. F. Cardone, M. G. Dainotti, S. Capozziello, and R. Willingale, Constraining cosmological parameters by gamma-ray burst X-ray afterglow light curves, MNRAS 408, 1181 (2010)
  * [22] M. G. Dainotti, V. F. Cardone, E. Piedipalumbo, and S. Capozziello, Slope evolution of GRB correlations and cosmology, MNRAS 436, 82 (2013)
  * [23] S. Postnikov, M. G. Dainotti, X. Hernandez, and S. Capozziello, Nonparametric Study of the Evolution of the Cosmological Equation of State with SNe Ia, BAO, and High-Redshift GRBs, ApJ 783, 126 (2014)
  * [24] M. G. Dainotti, M. Bogdan, A. Narendra, et al., Predicting the Redshift of $\gamma$-Ray-loud AGNs Using Supervised Machine Learning, ApJ 920, 118 (2021)
  * [25] M. G. Dainotti, V. Petrosian, M. Bogdan, et al., Gamma-ray Bursts as distance indicators through a machine learning approach, arXiv 2019, arXiv:1907.05074
  * [26] M. G. Dainotti, V. F. Cardone, and S. Capozziello, A time-luminosity correlation for $\gamma$-ray bursts in the X-rays, MNRASL 391, L79 (2008)
  * [27] M. G. Dainotti, R. Willingale, S. Capozziello, V. F. Cardone, and M. Ostrowski, Discovery of a Tight Correlation for Gamma-Ray Burst Afterglows with “Canonical” Light Curves, ApJL 722, L215 (2010)
  * [28] M. G. Dainotti, V. F. Cardone, S. Capozziello, M. Ostrowski, and R. Willingale, Study of possible systematics in the $L^{*}_{X}$ \- $T^{*}_{a}$ correlation of Gamma Ray Bursts, ApJ 730, 135 (2011)
  * [29] M. G. Dainotti, V. Petrosian, J. Singal, and M. Ostrowski, Determination of the intrinsic luminosity time correlation in the X-ray afterglows of gamma-ray bursts, ApJ 774, 157 (2013)
  * [30] M. G. Dainotti, R. Del Vecchio, N. Shigehiro, and S. Capozziello, Selection effects in Gamma-Ray Burst Correlations: Consequences on the Ratio Between Gamma-Ray Burst And Star Formation rates, ApJ 800, 31 (2015)
  * [31] M. G. Dainotti, V. Petrosian, R. Willingale, et al., Luminosity–time and luminosity–luminosity correlations for GRB prompt and afterglow plateau emissions, MNRAS 451, 3898 (2015)
  * [32] M. G. Dainotti, S. Postnikov, X. Hernandez, and M. Ostrowski, A fundamental Plane for Long Gamma-Ray Bursts with X-Ray plateaus, ApJL 825, L20 (2016)
  * [33] M. G. Dainotti, S. Nagataki, K. Maeda, S. Postnikov, and E. Pian, A study of gamma ray bursts with afterglow plateau phases associated with supernovae, A&A 600, A98 (2017)
  * [34] M. G. Dainotti, X. Hernandez, S. Postnikov, et al., A Study of the Gamma-Ray Burst Fundamental Plane, ApJ 848, 88 (2017)
  * [35] M. G. Dainotti, and L. Amati,Gamma-ray Burst Prompt Correlations: Selection and Instrumental Effects, PASP 130, 051001 (2018)
  * [36] M. G. Dainotti, R. Del Vecchio, R., and M. Tarnopolski, Gamma-Ray Burst Prompt Correlations, Adv. Astron. 2018, 4969503 (2018)
  * [37] M. G. Dainotti, A. L. Lenart, G. Sarracino, et al., The X-Ray Fundamental Plane of the Platinum Sample, the Kilonovae, and the SNe Ib/c Associated with GRBs, ApJ 904, 97 (2020)
  * [38] M. G. Dainotti, N. Omodei, G. P. Srinivasaragavan, et al., On the Existence of the Plateau Emission in High-energy Gamma-Ray Burst Light Curves Observed by Fermi-LAT, ApJS 255, 13 (2021)
  * [39] M. G. Dainotti, D. Levine, N. Fraija, and P. Chandra, Accounting for Selection Bias and Redshift Evolution in GRB Radio Afterglow Data, Galaxies 9, 95 (2021)
  * [40] R. Del Vecchio, M. G. Dainotti, and M. Ostrowski, Study of GRB Light-Curve Decay Indices in the Afterglow phase, ApJ 828, 36 (2016)
  * [41] M. G. Dainotti, V. Petrosian, and L. Bowden, Cosmological Evolution of the Formation Rate of Short Gamma-Ray Bursts with and without Extended Emission, ApJ 914, L40 (2021)
  * [42] S. Cao, M. G. Dainotti, and B. Ratra, Standardizing Platinum Dainotti-correlated gamma-ray bursts, and using them with standardized Amati-correlated gamma-ray bursts to constrain cosmological model parameters, arXiv 2022, arXiv:astro-ph.CO/2201.05245
  * [43] A. Rowlinson, B. P. Gompertz, M. G. Dainotti, et al., Constraining properties of GRB magnetar central engines using the observed plateau luminosity and duration correlation, MNRAS 443, 1779 (2014)
  * [44] N. Rea, M. Gullón, J. A. Pons, et al., Constraining The GRB-Magnetar Model by Means of the Galactic Pulsar Population, ApJ 813, 92 (2015)
  * [45] G. Stratta, M. G. Dainotti, S. Dall’Osso, X. Hernandez, and G. De Cesare, On the Magnetar Origin of the GRBs Presenting X-Ray Afterglow Plateaus, ApJ 869, 155 (2018)


${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\sigma}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>whether the sequence starting from
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\sigma}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>eventually drops to an empty stack, and if so, to which variable does it lead.
We can do this by, at each iteration and for each $n$, keeping track of the
set of configurations
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\omega}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>such that $\omega$ has size $n$ and the stack length has not dropped below $n$
since
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\omega}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>was last visited.

To conclude, we can decide whether a system is contractive by applying the
above procedure for every choice of
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\varepsilon}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>and
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\sigma}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>As there are only $|Q|(|\Delta|+1)$ (i.e., polynomially many) choices, the
total running time would be $\mathcal{O}(|Q|^{2}|\Delta|^{2})$. This can be
further reduced to $\mathcal{O}(|Q||\Delta|)$ since, after following the
procedure above for a given choice of
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\sigma}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>we can infer whether
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X^{\prime}}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\sigma^{\prime}}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>leads to an infinite sequence for all the
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X^{\prime}}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\sigma^{\prime}}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>that were visited, and thus exclude them from our subsequent analysis. In
other words, we can keep a lookup table that says, for every
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\sigma}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>whether it has been discovered before, and prune our search the next time we
revisit
${\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{X}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\langle}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color@cmyk@fill{1}{0.50}{0}{0}{\sigma}}{\color[rgb]{0,0.5,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,1}\pgfsys@color@cmyk@stroke{1}{0.50}{0}{0}\pgfsys@color<EMAIL_ADDRESS>∎


# Shape deformation of magnetically levitated fluid droplets

I. Sanskriti∗ D. Kim  J. Twamley<EMAIL_ADDRESS>Quantum Machines
Unit, Okinawa Institute of Science and Technology Graduate University, Onna,
Okinawa 904-0495, Japan

###### Abstract

Diamagnetic levitation can provide a completely passive method to support
materials against the pull of gravity, and researchers have levitated both
solids and fluids. Such levitation can be assisted by increasing the magnetic
susceptibility contrast by using a surrounding paramagnetic medium and through
buoyancy forces, known as magneto-Archimedean levitation. The magneto-
Archimedean levitation of solids has proved useful in chemistry and biology.
However, the levitation of fluid droplets has an additional interest because
the fluid droplet’s shape can deform. We perform experiments and simulations
to gauge the squashing or eccentricity of the static magnetically levitated
fluid droplet. By carefully characterizing all the parameters affecting the
droplet’s levitation, using image analysis to estimate the droplet’s
eccentricity, and using finite element adaptive simulations to find the lowest
energy droplet shape, we find good agreement between the simulations and
experimental results. As a potential application, we show that the droplet’s
eccentricity can be used to perform magnetic gradiometry with a potential
resolution of $S\sim 8\,{\rm nT/cm}$, over a volume of 10 mm3, which is
competitive with other room-temperature magnetic gradiometer techniques.

magnetic levitation, diamagnetic, liquid, shape

††preprint: AIP/123-QED

Studying the shape of liquid droplets has a long history. [1, 2] Developing a
deep understanding of the shape of liquid droplets and their dynamics can
reveal many properties, such as the influence of droplet shapes and dynamics
in radar sensing of rainfall activity. [3] The study of droplet shape is not
limited to fluid mechanics but also informs the study of processes over a wide
range of length scales, from the shapes of astrophysical bodies to the fission
of atomic nuclei.[4, 5] By levitating droplets, one can study their shape over
a long duration, and initial studies focused on the acoustic levitation of
droplets. [6, 7] Researchers even demonstrated the acoustic levitation of
liquid metal droplets. [8]

Diamagnetic levitation, which does not require active power, is highly
suitable for passive static levitation scenarios and has a long development
history. [9, 10, 11, 12] Researchers demonstrated the diamagnetic levitation
of various organic fluids [13], but this required a bitter solenoid or
superconducting magnet with the approximate magnetic field of 15-20 T.
Researchers could also magnetically levitate liquid helium droplets, but it
required $|B(dB/dz)|$ of $20.7$ T${}^{2}/$cm. [14] Later, by combining
buoyancy forces and a paramagnetic medium, researchers levitated fluids at
even lower magnetic field strengths of 10 T and named the method magneto-
Archimedes levitation.[15] These studies opened up a new dimension of
diamagnetic levitation and its potential applications in various fields,
including chemistry, biology, materials science, and biochemistry. [16] The
development of magneto-Archimedean levitation of solids and liquids led
researchers to study the dynamics, shape, and size of levitated static drops
or levitated spinning drops affected by magnetic forces. [17, 18, 19, 20, 21,
22] Researchers also explored this technique to separate biological materials,
[23] measure densities of solids and water immiscible organic liquids, [24,
25] develop density-based sensors to characterize food and water, [26] and
separate and identify illicit drugs. [27]

Many of these works focused on studying the shape dynamics of the fluid
droplet, while others focused on the levitation height of the fluid droplet.
In this work, we focus on determining the deformation in the shape of the
magneto-Archimedean levitated droplet due to the combination of magnetic,
buoyancy, gravity, and surface tension forces. There have been only a few
previous studies on the magnetic shape deformation of diamagnetically fluid
drops, but these either consider pendant drops, [18, 28] or air bubble
magnetically levitating within liquid oxygen in a large superconducting
electromagnet. [29]

In this work, we explore if the magnetic forces in a magneto-Archimedean fluid
levitation setup can be engineered to affect a significant deformation of the
static shape of the levitated droplet. Instead of superconducting magnets, we
have used permanent rare-earth magnets and a diamagnetic organic fluid droplet
levitated in a paramagnetic aqueous salt medium. Here, we demonstrate that the
levitated droplet exhibits significant squashing in its shape and becomes
oblate by adjusting the position of the trapping magnets due to the increased
strength of the magnetic forces on the droplet in an anisotropic fashion. We
performed experiments to measure the shape deformations of the droplet as we
altered the magnetic forces, and these agreed with numerical simulation.

Our findings indicated that the droplet eccentricity depends on the magnetic
field gradients at the droplet’s location, which led our investigation to
whether the eccentricity can be used as a sensor for magnetic gradients, e.g.,
work as a magnetic gradiometer. There exist already various types of magnetic
field sensors or gradiometers, each with its strengths and weaknesses. Through
simulations, we found that the estimation of the droplet eccentricity can
yield a magnetic gradiometer whose sensitivity is competitive with many
existing magnetic gradiometer technologies, operates at room temperature, and
is compact and robust.

In the following, we describe the physics of the magneto-Archimedean droplet
trap, outline the experimental setup, and then work to carefully quantify each
of the parameters and variables that define the trapping forces. As a check,
we examined theoretically and experimentally the equilibrium droplet
levitation heights as we varied the magnetic forces (by moving the magnets),
which agreed well with previously published works. [24, 25] We performed
experimental measurements and numerical simulations of the eccentricity of the
static droplet’s shape. Finally, we examined the potential for the shape
deformation to act as a magnetic gradient sensor.

For describing the physics of trapping, we considered an infinitesimal element
of the diamagnetic fluid of volume $V$, located at position $\vec{r}$. We can
express the total net force density it experiences as [24],
$\vec{F}(\vec{r})/V\equiv\vec{f}=-\vec{\nabla}u(\vec{r})$, and the potential
density $u(\vec{r})$ is given as

$u(\vec{r})=(\rho_{s}-\rho_{m})g\,z-\frac{1}{2\mu_{0}}(\chi_{s}-\chi_{m})|\vec{B}(\vec{r})|^{2}\;\;,$
(1)

where $(\rho_{s},\chi_{s})$, is the (density [${\rm kg/m^{3}}$], volume
magnetic susceptibility [unit-less]), of the sample (also called the droplet
fluid), while $(\rho_{m},\chi_{m})$ is that of the surrounding medium (or
paramagnetic fluid). The diamagnetic sample fluid, 3-chlorotoluene, has known
values of $\rho_{s}$ and $\chi_{s}$ and was used as received from the
suppliers. For the paramagnetic medium, $\rm{MnCl_{2}\cdot 4H_{2}O}$, the
density and magnetic susceptibility will vary depending on the prepared molar
concentration of dissolved paramagnetic salts in an aqueous solution. The
magnetic field $\vec{B}(\vec{r})$ [T] is generated by two opposing ring
magnets whose inter-magnet separation is $d$ [m], and $z$ [m] is the
$z-$component of $\vec{r}$, the height of the infinitesimal volume element
above the top surface of the bottom magnet. We take the origin of the axes as
the axial center of the top surface of the bottom ring magnet.

We can estimate the equilibrium height, $h$, of the levitated droplet by
solving

$f_{z}(\vec{r}=(0,0,z=h))=0,\;\;\Rightarrow\;\;\left.\frac{d}{dz}\,u(\vec{r}=\\{0,0,z\\})\right|_{z=h}=0\;\;.$
(2)

This approximation is justified when the droplet volume is infinitesimal so
that $\int_{V}f_{z}(\vec{r})\approx 0$. In our case, the droplet size is
comparable to the inter-magnetic separation $d$, and we find that the magnetic
potential around the droplet is non-Gaussian. Therefore, we can use the
following to obtain a more precise estimation of the equilibrium droplet
height, $h$, by solving

$F_{z}(\vec{r}=(0,0,z=h))=0,\;\;\Rightarrow\;\;\left.\frac{d}{dz}\left(\,\int_{V(z,r)}u(\vec{r})\,dV\right)\right|_{z=h}=0\;\;.$
(3)

where $V(z,r)$ denotes a spherical integration volume of radius $r$, the
droplet radius, where the sphere is centered at the position
$\vec{r}=\\{0,0,z\\}$. As an example, we show the potential density
$u(r,z^{*})$, in Fig. 1, where we set $z^{*}=z-h$.

To evaluate the equilibrium shape of the droplet when it is levitating in the
potential density $u(\vec{r})$, we also require the value of the interfacial
tension (IFT), $\gamma$ [N/m], between the droplet fluid and the medium fluid.
We see that the small changes in $\gamma\sim\pm 2$ mN/m result in changes in
the droplets’ eccentricity by a few parts in a hundred. However, if the change
in $\gamma\sim\pm 30$ mN/m, the change in eccentricity can be a few parts in
ten.

Figure 1: Example potential density plot. We graph the Archimedean Magnetic
potential density given in Eq. 1, in axial coordinates, which includes the
effects of the magnetic, gravitational, and buoyancy forces. We observe a
minimum that traps the diamagnetic droplet in space. Here the magnet
separation is set as $d=7\,{\rm mm}$, with the droplet volume
$V_{drop}=10\,{\rm\mu L}$, and radius $r\sim 1.33\,{\rm mm}$, with a $3\,{\rm
M}$ concentration of the paramagnetic medium. Solving Eq. 2 yields a
levitation height of $h\sim 5.25\,{\rm mm}$, and we choose $z^{*}=z-h$ to
center the plot at this height. The contours are spaced by $\Delta
u=0.5\,{\rm[J/m^{3}]}$. The red disc indicates the approximate position and
dimensions of the diamagnetic droplet. Although the potential exhibits
significant oblateness, the deformation of the droplet towards these contours
is resisted by the high interfacial tension. In this example, we have chosen:
magnet dimensions [OD, ID, H]: [38.6, 19.1, 30.9] ${\rm mm}$, the
magnetization of the magnets are assumed to be identical $M\sim 1.32\,{\rm
T}$. The droplet fluid is considered to be pure 3-chlorotoluene, and the
medium fluid is a $3\,{\rm M}$ concentration of aqueous $\rm MnCl_{2}\cdot
4H_{2}O$, with the associated physical properties as listed in [Table S1, SI].

We consider the setup depicted in Fig. 2, which is similar to that in [25],
consisting of two opposed NdFeB (N40 - Neomag) ring magnets with dimensions
$38.6\,{\rm mm}$ (outer diameter) $\times\,19.1\,{\rm mm}$ (inner diameter)
$\times\,30.9\,{\rm mm}$ (height), surrounding a rectangular quartz cuvette
containing a paramagnetic solution of $\rm MnCl_{2}\cdot 4H_{2}O$, prepared
with various concentrations. A fixed-volume droplet of diamagnetic liquid
(3-chlorotoluene) was inserted into the paramagnetic medium and was trapped by
a combination of buoyancy and magnetic forces. As indicated by the green arrow
in Fig. 2(c), we can adjust the position of the upper magnet and, in
particular, the inter-magnet separation $d$ via a robust micrometer stage.
This stage has to withstand considerable repulsive forces generated between
the magnets when $d$ is reduced. The setup shown also has affixed a digital
micrometer to record the inter-magnet separations accurately. For measuring
the levitation height and capturing the image of the levitated droplet between
the two magnets, a microscope (Dino-lite Edge 3.0), along with lights, was
mounted to illuminate the levitated droplets.

Figure 2: Schematic of magnetic levitation and squashing of droplet: (a) a
small (1-2 mm diameter) droplet of a diamagnetic liquid (round sphere) is
trapped by a combination of magnetic and buoyancy forces within a rectangular
glass cuvette, containing a paramagnetic liquid (clear). The cuvette sits
between two opposing NdFeB ring magnets, producing an anti-Helmholtz field.
The separation distance, $d$, between these magnets can be varied. (b) detail
of the squashed levitated droplet (green - exaggerated for clarity), with a
levitation height $h$, magnet separation $d$, and major (minor) droplet radii
$a\,(b)$, and eccentricity $e$; (c) photo of the experimental setup
illustrating the ring magnets and reinforced micrometer positioning
arrangement of the upper ring magnet indicated by a green arrow, and
micrometer gauge to measure the magnet separation. Not shown are the rigs to
hold the microscope camera to measure the droplet shape, lights to illuminate
the droplet, and micrometer Hall-probe scanning stage to measure the trapping
$B$ fields.

To precisely utilize Eq. 1, we must have precise information about the
densities and magnetic susceptibilities of the sample and medium and the
magnetic fields. The preparation of various concentrations of the paramagnetic
medium solution is described in [Sec III, SI]. The density of the medium
solution $\rho_{m}(T,c)$ is a function of temperature and concentration of the
paramagnetic salt $c$. We measured $\rho_{m}$ for various values of $c$ using
a pycnometer [Sec IV, SI]. The experimental values of density of (1 M, 2 M,
and 3 M) aqueous solutions of MnCl${}_{2}\cdot$4H2O at 25 ∘C were measured to
be (1.097 $\pm$ 0.004, 1.195 $\pm$ 0.005, 1.293 $\pm$ 0.005) g/mL respectively
that compared well with the empirical formula [Eq S4, SI]. The magnetic
susceptibility of the 3-chlorotoluene $\chi_{s}$ is taken from the literature,
while $\chi_{m}(T,c)$ for the paramagnetic medium is a function of temperature
and salt concentration. Values for $\chi_{m}(T,c)$ are calculated using the
Curie-Weiss law [Sec V, SI]. The potential density is also a function of the
magnetic fields generated by the two ring magnets. By setting the ring magnets
to have a separation of $d=12$ mm, we experimentally measured the $B_{z}(z)$
field along the $z-$axis between the two ring magnets using a Hall probe
magnetometer. By fitting the experimental data to theory, we estimated the
magnetization of the top and bottom ring magnets ($M_{top}$ = 1.36 T and
$M_{bottom}$ = 1.3 T), which are quite comparable to those quoted by the
manufacturer [Sec VI, SI]. These values were used for all subsequent
simulations and modeling. To check the accuracy of the parameters and setup,
we measured the levitation height of the droplet for various inter-magnet
separations $d\sim 6-12$ mm and compared the experimental data with
simulations and found quite good agreement [Sec VII, Fig S5, SI]. The slight
systematic upwards shift of the experimental data compared to the theoretical
predictions may be due to a slight tilting of the movable ring magnet away
from an exactly horizontal orientation.

The droplet’s shape will be greatly affected by the value of the interfacial
tension (IFT) $\gamma$ between the droplet fluid and the surrounding medium.
Hence, we determined the IFT between 3-chlorotoluene and a 3 M aqueous
solution of MnCl${}_{2}\cdot$4H2O, $\gamma_{CT/MnCl_{2}}$, experimentally with
the help of Pendant drop Tensiometry. [30] The value for
$\gamma_{CT/MnCl_{2}}$ was found to be $45.23\pm 0.55\,{\rm mN/m}$. This value
is used later in simulations of the droplet shape. The experimental details,
along with the figure of the pendant droplet, are given in [Sec VIII, SI].

Figure 3:  Image analysis of a droplet shape for inter-magnet separation $d=9$
mm. (a) [Left] we schematically depict how the radius of the droplet depends
on the angle $\theta$, (a) [Right] actual microscope image of the levitated
drop, (b) the analysis of the droplet image in (a) [Right], where the blue
points label points on the edge of the droplet and the data is repeated
outside the fundamental domain $\theta\in[0,2\pi]$. The red curve is a fitted
function for $r(\theta)$, using a high-order Fourier series. From the minimum
and maximum of $r(\theta)$ and the spread of values around these, we estimate
the eccentricity of this droplet to be $e=0.093\pm 0.018$.

We performed experiments and image analysis using Mathematica to determine the
eccentricity of the droplet with varying inter-magnet separation $d$. A
digital microscope was used to take close-up photos of the levitated droplet.
The lighting conditions were adjusted to capture photos where the droplet edge
exhibited high contrast with its background. Two side light sources
illuminated the levitated droplet. To achieve high droplet edge contrast in
the images, we placed a black square piece of hardboard with a circular
central white spot (2 cm in diameter), which is placed in the background
behind the droplet as seen by the microscope. This leads to an image where the
background and the droplet interior are white, but the droplet edge is dark.
This results in a high-contrast image of the edge of the droplet. As the
inter-magnet separation distance was changed, several microscope images of the
droplet were taken. Image processing was performed on the high-contrast edge
of the droplet to identify the coordinates of points. From this, an estimation
of the eccentricity of the droplet image was obtained. The droplet image was
processed by first binarizing and blackening the interior. This created a
higher contrast for further processing for edge detection. Then, the 2D
coordinates of points on the edge were isolated. The center point of the
droplet was also identified during image processing. Next, the radial distance
from this center point to each point $i$ on edge and angle to the horizontal
was found, e.g., $(r_{i},\theta_{i})$, where $\theta_{i}\in[-\pi,+\pi]$. The
experimental data was found to be periodic in $\theta$ with period $2\pi$, and
we used this periodicity to expand the domain of $\theta\in[-3\pi,+3\pi]$. We
dropped some data points clustered around $\theta\sim m\pi/2$, where $m$ is an
integer, as the finite resolution of the camera image, which is in a cartesian
grid, does not map very smoothly under the polar transformation at these
angles [28].We fit a Fourier expansion to this periodic data as
$r(\theta)=\sum_{k=0}^{30}\,[A_{k}\sin(k\pi\theta+\phi)+B_{k}\cos(k\pi\theta+\phi)]$,
fitting for the amplitudes $A_{k},\,B_{k}$, and phase shift $\phi$. The latter
parameter describes any slight mismatch of the microscope photo’s vertical
axis from the true lab vertical axis. From this fitting, we found the major
axis $a=\max(r(\theta))$ and minor axis $b=\min(r(\theta))$ of the droplet. We
estimated the error in these quantities, $(\delta a,\delta b)$, by the
standard deviation of $r_{i}$ around these locations. An example plot of
$r_{i}$ vs. $\theta_{i}$, for a drop with $d=9$ mm, is shown in Fig. 3. From
this, we observe that these errors are not small compared to the mean value of
the eccentricity. One can find the error in the derived eccentricity
$e=\sqrt{1-(b/a)^{2}}$, as

$\delta e=\left(\frac{b^{2}}{a^{3}}\right)\frac{\delta
a}{e}+\left(\frac{b}{a^{2}}\right)\frac{\delta b}{e}\;\;.$ (4)

As the inter-magnet separation distance $d$, was varied from $d=7,8,9,10,11$
mm, for each separation, we took several photo-micro-graphs of the droplet,
performed the above-described image analysis, and for each image, obtained an
estimate for the eccentricity, and its error $e\pm\delta e$. In the cases when
$d>11$ mm and $d<7$ mm, the image analysis failed as the droplet was partially
obscured by the magnets [Sec IX, SI].

In the subsequent sections, we describe how we model and predict the shape of
the levitated droplet numerically, in particular, estimate the droplet’s
eccentricity. Since the magnetic fields produced by the ring magnets are
axially symmetric, the droplet shape will also be axially symmetric (i.e., the
shape will be independent of $\phi$), and we can thus work in a reduced two-
dimensional $x-z$ plane. We numerically sampled $u(\vec{r})$ in a 2D region
surrounding the droplet on this plane using Eq. 1, where the $B$ field
produced by the magnets is simulated using MagPyLib. [31]The shape of the
droplet is that which minimizes the surface energy in the presence of
$u(\vec{r})$. Two limiting cases are when a) if the interfacial tension
vanishes $\gamma=0$, the shape of the droplet exactly follows the
equipotentials of the potential density $u(\vec{r})$, or b) if
$\gamma\sim\infty$, then the shape of the droplet is a sphere. However, in the
case of finite $\gamma$, estimating the droplet’s shape is not
straightforward. Such surface energy minimization calculations can be
performed using the Surface Evolver (SE). [32, 33, 34] SE can estimate droplet
shape through iteration and refinement, including additional force effects and
a spatially varying potential energy density. It also considers the
interfacial tension and the volume of the droplet enveloped by the surface.

The Surface Evolver program incorporates spatial potentials as analytic
formulae. Thus we sampled the potential density $u(\vec{r})$ around the small
droplet and fitted as a polynomial function of degree 2 in the coordinates $x$
and $z$, where the origin is taken as the droplet’s equilibrium levitation
height.We then called the SE program using this potential, incorporating the
interfacial surface tension and droplet volume used in all experiments. After
appropriate iterations and refinements, vertices of the evolved shape are
exported for further analysis. From the SE exported data, we calculated the
radii of all surface points to the geometric center of the 3D surface and set
the major and minor radii of the oblate droplet to be the maximum and minimum
radii found. From this, we computed the droplet eccentricity. As a study, we
performed simulations for varying magnet separation distance
$d\in\\{4,16\\}\,{\rm mm}$ and for various test values of interfacial
tensions, and these results are reported in [Fig S10, SI]. As one might
expect, the droplet’s eccentricity is increased with lower interfacial
tensions $\gamma$ and lower magnet separation distance $d$. We also observed
that the simulated eccentricity value tends to be a constant for large
separations $d$ greater than 10 mm. In this case, the magnetic forces from a
single-ring magnet are sufficient to trap the droplet. Since in magneto-
statics, one must have
$\vec{\nabla}\cdot\vec{B}(\vec{r})=0,\,\forall\,{\vec{r}}$, these trapping
forces cannot be isotropic, and there is some residual squashing of the
droplet. We simulated the droplet’s shape in the actual experiments we
performed. We plotted $e(d)$, for $V=$ 10 ${\rm\mu L}$, using the
experimentally determined value of the interfacial tension $\gamma=45\,{\rm
mN/m}$, as the black curve in Fig. 4. We also developed extensive methods to
estimate the accuracy of the numerical simulations as SE’s algorithm is
somewhat stochastic [Sec X, SI]. We found that SE returns an estimate for the
eccentricity, which is always lower than the actual value. For the experiment,
we depict the errors in the numerical simulation in Fig. 4 as a grey region
below the black solid line.

Figure 4: Comparison between simulation and experimental data of eccentricity.
Parameters used for the simulation: $\gamma$ = 45 mN/m, $V$=10 $\mu$L, $c$ = 3
M. We plot the mean eccentricity (plus markers) and errors for the
experimentally determined droplet eccentricity for inter-magnet separations
$d=7,8,9,10,11$ mm. The numerical simulation is the solid black line which is
an upper bound for the estimated numerical eccentricity. The associated
numerical error with this computation is shaded in gray. The simulations and
experimental data agree except for the purple outlier at $d=10$ mm, which
suffered from poor background non-uniformity.

In Fig. 4, we plot both the experimental and numerical simulation results for
the droplet eccentricities, and this summarizes the main result of our
inquiry. We note that all the droplets are oblate, even for inter-magnet
separations where the droplet is trapped by only one magnet. We observed that
the numerical simulations, including error, agree with the experimental
measurements predicting a slight increase in eccentricity as the inter-magnet
separation is decreased below 9 mm, which is apparent in the experimental data
for $d>7$ mm [See Sec XI, SI for details].

We also investigated the potential of utilizing eccentricity as a sensor of
the local magnetic field gradient. Different magnetic gradiometry methods have
varying sensitivities and applications due to their characteristics. The Hall
sensor has relatively low sensitivity, but due to its versatility, it is
commonly used in daily life, such as in smartphones. On the other hand,
sensors using microelectromechanical systems (MEMS), superconductors, and
atoms have high sensitivities, but their usage is limited due to the technical
specifics required by these techniques. From Eq. 1, we observe that the
magnetic squashing force, and thus the eccentricity of the droplet, depends on
the gradient of the local magnetic fields. We ask whether the eccentricity can
be used as a sensor of the magnetic field gradient $S$ [T/m], e.g. if we
measure the major and minor radii of the droplet to a preset precision, e.g.,
$\delta r~{}\sim 1$ pm, what is the minimum resolvable change in the magnetic
field gradient that can be sensed by estimating the eccentricity? To estimate
$S$, we make use of the SE and numerically estimate

$S(d_{0})=\left.\frac{\text{d}B_{z}^{\prime}(d)}{\text{d}a(d)}\right|_{d=d_{0}}\times
a_{min}\;\;.$ (5)

where we assume the resolution to measure $a$ is $a_{min}\sim 10^{-12}$ m,
which is possible using laser interferometry (for example, the Picoscale
interferometer). Using an interfacial tension value of 1 mN/m simulations
indicated one can achieve $S\sim$ 8 nT/cm [Sec XII, SI]. The predicted values
of $S$ are comparable with other sensitive magnetic gradiometer technologies,
which can operate in ambient conditions and have a volume of 10 mm3.

The authors acknowledge support by the Okinawa Institute of Science and
Technology Graduate University. The authors also acknowledge technical
assistance from P. Kennedy from the OIST Engineering Section.

The data and simulation codes that support the findings of this study are
available from the corresponding author upon reasonable request.

## References

  * Rayleigh [1879] L. Rayleigh, “VI. On the capillary phenomena of jets,” Proceedings of the Royal Society of London 29, 71–97 (1879).
  * Lamb [1892] H. Lamb,  _Hydrodynamics_, 1st ed. (Cambridge University Press, Cambridge [Eng.] :, 1892).
  * Hill and Eaves [2012] R. J. A. Hill and L. Eaves, “Shape oscillations of an electrically charged diamagnetically levitated droplet,” Applied Physics Letters 100, 114106 (2012).
  * Holgate and Coppins [2018] J. T. Holgate and M. Coppins, “Shapes, stability, and hysteresis of rotating and charged axisymmetric drops in a vacuum,” Physics of Fluids 30, 064107 (2018).
  * Butler [2020] S. L. Butler, “Equilibrium shapes of two-phase rotating fluid drops with surface tension,” Physics of Fluids 32, 012115 (2020).
  * Lee, Anilkumar, and Wang [1991] C. P. Lee, A. V. Anilkumar, and T. G. Wang, “Static shape and instability of an acoustically levitated liquid drop,” Physics of Fluids A: Fluid Dynamics 3, 2497–2515 (1991).
  * Yarin, Pfaffenlehner, and Tropea [1998] A. L. Yarin, M. Pfaffenlehner, and C. Tropea, “On the acoustic levitation of droplets,” Journal of Fluid Mechanics 356, 65–91 (1998).
  * Xie  _et al._ [2002] W. J. Xie, C. D. Cao, Y. J. Lü, and B. Wei, “Levitation of Iridium and Liquid Mercury by Ultrasound,” Physical Review Letters 89, 104304 (2002).
  * Thomson [1850] W. Thomson, “Remarks on the forces experienced by inductively magnetized ferromagnetic or diamagnetic non-crystalline substances,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 37, 241–253 (1850).
  * Braunbek [1939] W. Braunbek, “Freies Schweben diamagnetischer Körper im Magnetfeld,” Zeitschrift für Physik 112, 764–769 (1939).
  * Berry and Geim [1997] M. V. Berry and A. K. Geim, “Of flying frogs and levitrons,” European Journal of Physics 18, 307–313 (1997).
  * Simon and Geim [2000] M. D. Simon and A. K. Geim, “Diamagnetic levitation: Flying frogs and floating magnets (invited),” Journal of Applied Physics 87, 6200–6204 (2000).
  * Beangnon and Tournier [1991] E. Beangnon and R. Tournier, “Levitation of organic materials,” Nature 349, 470–470 (1991).
  * Weilert  _et al._ [1996] M. A. Weilert, D. L. Whitaker, H. J. Maris, and G. M. Seidel, “Magnetic levitation and noncoalescence of liquid helium,” Physical Review Letters 77, 4840–4843 (1996).
  * Ikezoe  _et al._ [1998] Y. Ikezoe, N. Hirota, J. Nakagawa, and K. Kitazawa, “Making water levitate,” Nature 393, 749–750 (1998).
  * Ge  _et al._ [2020] S. Ge, A. Nemiroski, K. A. Mirica, C. R. Mace, J. W. Hennek, A. A. Kumar, and G. M. Whitesides, “Magnetic Levitation in Chemistry, Materials Science, and Biochemistry,” Angewandte Chemie International Edition 59, 17810–17855 (2020).
  * Whitaker  _et al._ [1998] D. L. Whitaker, C. Kim, C. L. Vicente, M. A. Weilert, H. J. Maris, and G. M. Seidel, “Shape Oscillations in Levitated He II Drops,” Journal of Low Temperature Physics 113, 491–499 (1998).
  * Katsuki  _et al._ [2007] A. Katsuki, K. Kaji, M. Sueda, and Y. Tanimoto, “Magnetic Field Effects on Shape and Size of Pendent Water Drop,” Chemistry Letters 36, 306–307 (2007).
  * Hill and Eaves [2008] R. J. A. Hill and L. Eaves, “Nonaxisymmetric Shapes of a Magnetically Levitated and Spinning Water Droplet,” Physical Review Letters 101, 234501 (2008).
  * Hill and Eaves [2010] R. J. A. Hill and L. Eaves, “Vibrations of a diamagnetically levitated water droplet,” Physical Review E 81, 056312 (2010).
  * Temperton, Hill, and Sharp [2014] R. H. Temperton, R. J. A. Hill, and J. S. Sharp, “Mechanical vibrations of magnetically levitated viscoelastic droplets,” Soft Matter 10, 5375–5379 (2014).
  * Liao and Hill [2017] L. Liao and R. J. A. Hill, “Shapes and Fissility of Highly Charged and Rapidly Rotating Levitated Liquid Drops,” Physical Review Letters 119, 114501 (2017).
  * Hirota  _et al._ [2004] N. Hirota, M. Kurashige, M. Iwasaka, M. Ikehata, H. Uetake, T. Takayama, H. Nakamura, Y. Ikezoe, S. Ueno, and K. Kitazawa, “Magneto-Archimedes separation and its application to the separation of biological materials,” Physica B: Condensed Matter 346-347, 267–271 (2004).
  * Mirica  _et al._ [2009] K. A. Mirica, S. S. Shevkoplyas, S. T. Phillips, M. Gupta, and G. M. Whitesides, “Measuring Densities of Solids and Liquids Using Magnetic Levitation: Fundamentals,” Journal of the American Chemical Society 131, 10049–10058 (2009).
  * Ge and Whitesides [2018] S. Ge and G. M. Whitesides, “"Axial" Magnetic Levitation Using Ring Magnets Enables Simple Density-Based Analysis, Separation, and Manipulation,” Analytical Chemistry 90, 12239–12245 (2018).
  * Mirica  _et al._ [2010] K. A. Mirica, S. T. Phillips, C. R. Mace, and G. M. Whitesides, “Magnetic Levitation in the Analysis of Foods and Water,” Journal of Agricultural and Food Chemistry 58, 6565–6569 (2010).
  * Abrahamsson  _et al._ [2020] C. K. Abrahamsson, A. Nagarkar, M. J. Fink, D. J. Preston, S. Ge, J. S. Bozenko, and G. M. Whitesides, “Analysis of Powders Containing Illicit Drugs Using Magnetic Levitation,” Angewandte Chemie International Edition 59, 874–881 (2020).
  * Poulose  _et al._ [2022] S. Poulose, J. A. Quirke, P. Stamenov, M. E. Möbius, and J. M. D. Coey, “Deformation and necking of liquid droplets in a magnetic field,” Physics of Fluids 34, 112116 (2022).
  * Duplat and Mailfert [2013] J. Duplat and A. Mailfert, “On the bubble shape in a magnetically compensated gravity environment,” Journal of Fluid Mechanics 716, R11 (2013).
  * Berry  _et al._ [2015] J. D. Berry, M. J. Neeson, R. R. Dagastine, D. Y. Chan, and R. F. Tabor, “Measurement of surface and interfacial tension using pendant drop tensiometry,” Journal of Colloid and Interface Science 454, 226–237 (2015).
  * Ortner and Coliado Bandeira [2020] M. Ortner and L. G. Coliado Bandeira, “Magpylib: A free Python package for magnetic field computation,” SoftwareX 11, 100466 (2020).
  * Brakke [1992] K. A. Brakke, “The Surface Evolver,” Experimental Mathematics 1, 141–165 (1992).
  * Brakke and Sullivan [1997] K. A. Brakke and J. M. Sullivan, “Using Symmetry Features of the Surface Evolver to Study Foams,” in  _Visualization and Mathematics_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1997) pp. 95–117.
  * Zhang  _et al._ [2016] Y. Zhang, D. Chatain, S. L. Anna, and S. Garoff, “Stability of a compound sessile drop at the axisymmetric configuration,” Journal of Colloid and Interface Science 462, 88–99 (2016).