# New small gaps between squarefree numbers

Tsz Ho Chan

###### Abstract

In this paper, we show that, for some constant $C>0$, the interval
$(x,x+Cx^{5/26}]$ always contains a squarefree number when $x$ is sufficiently
large (in terms of $C$). Our improvement comes from establishing asymptotic
relations between the shifts $a$ and $b$ when $mn^{2}\approx(m-a)(n+b)^{2}$ We
apply them to study quadruples $(m+a_{1})(n-b_{1})^{2}\approx
mn^{2}\approx(m-a_{2})(n+b_{2})^{2}\approx(m-a_{2}-a_{3})(n+b_{2}+b_{3})^{2}$
and generalize Roth differencing and Filaseta-Trifonov differencing by
allowing $b_{1}$ to be different from $b_{3}$. We also introduce a new
differencing and exploit the interplay among these three differencings.

##  1 Introduction and Main Result

A positive integer is squarefree if it is not divisible by the square of any
prime numbers. For example, $6$ is squarefree but $12$ is not. It is well-
known that

$\mathop{\sum_{n\leq x}}_{n\text{
squarefree}}1=\frac{6}{\pi^{2}}x+O(x^{1/2}).$ (1)

Thus, the set of squarefree numbers has positive density and it is believed
that gaps between successive squarefree numbers is small.

###### Conjecture 1

For any $\epsilon>0$ and $x$ sufficiently large (in terms of $\epsilon$), the
short interval $(x,x+x^{\epsilon}]$ always contains a squarefree number.

Various authors have studied this question. Firstly, Fogel [6] showed that the
conjecture is true for $\epsilon>2/5$ by studying average of arithmetic
functions over short intervals. This was greatly improved by Roth [10] to
$\epsilon>1/4$ through an elementary argument and $\epsilon>3/13$ by
incorporating van der Corput’s method. Later, Richert [9], Rankin [8], Schmidt
[11], Graham and Kolesnik [7], Filaseta [1, 2], Trifornov [12, 13] and
Filaseta and Trifonov [3] improved the result by exponential sum and
elementary techniques. The current best result was by Filaseta and Trifonov
[4]:

###### Theorem 1

There exists a constant $C_{0}>0$ such that, for $x$ sufficiently large, the
interval $(x,x+C_{0}x^{1/5}\log x]$ contains a squarefree number.

In this paper, we are going to improve Theorem 1 to

###### Theorem 2

There exists a constant $C_{0}>0$ such that, for $x$ sufficiently large, the
interval $(x,x+C_{0}x^{5/26}]$ contains a squarefree number.

Acknowledgement The author would like to thank Prof. Hugh Montgomery and Prof.
Trevor Wooley for introducing him to this subject as a graduate student. He
would also like to thank Prof. Michael Filaseta and Prof. Ognian Trifonov for
encouragement and suggestions as well as pointing out mistakes in the author’s
first attempt.

##  2 Some Preparations

Throughout this paper,

$\displaystyle c\text{'s are sufficiently small positive constants}.$
$\displaystyle C\text{'s are sufficiently large positive constants}.$
$\displaystyle L>0\text{ is a large parameter}.$ $\displaystyle\lambda=1+1/L.$
$\displaystyle x\text{ is assumed to be sufficiently large in terms of
}c\text{'s and }C\text{'s}.$ $\displaystyle p\text{ denotes a prime number}.$
$\displaystyle|S|\text{ denotes the cardinality of a set }S.$
$\displaystyle\|\theta\|\text{ denotes the distance between $\theta$ and its
nearest integer}.$ $\displaystyle\theta=5/26.$ $\displaystyle\delta=0.0751.$
$\displaystyle h=C_{0}x^{\theta}.$ $\displaystyle f(x)\ll g(x)\text{ or
}g(x)\gg f(x)\text{ means that there is a constant }C\text{ independent of
}C_{0}\text{ and }L.$ $\displaystyle\text{ such that }|f(x)|\leq Cg(x).$
$\displaystyle f(x)\asymp g(x)\text{ means that }f(x)\ll g(x)\text{ and
}g(x)\ll f(x).$ $\displaystyle f(x)=o(g(x))\text{ means that
}\lim_{x\rightarrow\infty}f(x)/g(x)=0.$ $\displaystyle mn^{2}\text{ is some
number in the interval }(x,x+h].$ $\displaystyle n\text{ is typically in the
interval }(N,\lambda N].$ $\displaystyle m\text{ typically has size
}M=x/N^{2}.$ $\displaystyle R=c_{0}N^{4/3}/x^{1/3}.$ $\displaystyle
Q=\max(R,x^{1/5-\delta}).$

We want to show that the interval $(x,x+h]$ contains a squarefree number. We
may assume that $h\leq 0.1x$ for otherwise the result follows from (1). Our
initial treatment follows that of Filaseta and Trifonov [4]. Let $S$ denote
the set of integers in $(x,x+h]$ which are not squarefree. Then any member of
$S$ must be of the form $mn^{2}$ for some integers $m\geq 1$ and $n>1$. Let
$M_{n}$ denote the number of multiples of $n^{2}$ in the interval $(x,x+h]$.
Then, for $x$ sufficiently large, we have $M_{n}=0$ if $n>2\sqrt{x}$. Since
every number greater than one has a prime factor,

$\displaystyle|S|$ $\displaystyle\leq\sum_{p\leq 2\sqrt{x}}M_{p}=\sum_{p\leq
h}M_{p}+\sum_{h<p\leq ch^{5/2}/x^{1/4}}M_{p}+\sum_{ch^{5/2}/x^{1/4}<p\leq
x^{2/3}\log^{2}x/h^{5/3}}M_{p}+\sum_{x^{2/3}\log^{2}x/h^{5/3}<p\leq
2\sqrt{x}}M_{p}$ $\displaystyle=:S_{1}+S_{2}+S_{3}+S_{4}$ (2)

where $c>0$ is some small constant and $C^{\prime}>0$ is some large constant.
Clearly, $M_{p}\leq h/p^{2}+1$ and

$S_{1}\leq\sum_{p\leq
h}\Bigl{(}\frac{h}{p^{2}}+1\Bigr{)}<h\sum_{n=2}^{\infty}\frac{1}{n^{2}}+\pi(h)=h\Bigl{(}\frac{\pi^{2}}{6}-1\Bigr{)}+\pi(h)\leq\frac{2}{3}h$

by a Chebyshev estimate. Define

$S(t_{1},t_{2}):=\\{n\in(t_{1},t_{2}]:\text{ there exists an integer $m$ such
that }mn^{2}\in(x,x+h]\\}$

Suppose $n\in(h,2\sqrt{x}]$. Then there can be at most one multiple of $n^{2}$
in the interval $(x,x+h]$ as $n^{2}>h$. Hence, for $n\in(h,2\sqrt{x}]$,
$M_{n}=1$ if and only if $n\in S(h,2\sqrt{x})$ and $M_{n}=0$ if and only if
$n\not\in S(h,2\sqrt{x})$. Therefore,

$S_{2}\leq|S(h,ch^{5/2}/x^{1/4})|,\;\;S_{3}\leq|S(ch^{5/2}/x^{1/4},x^{2/3}\log^{2}x/h^{5/3})|,\;\;S_{4}\leq|S(x^{2/3}\log^{2}x/h^{5/3},2\sqrt{x})|,$

and

$|S|\leq\frac{2}{3}h+|S(h,ch^{5/2}/x^{1/4})|+|S(ch^{5/2}/x^{1/4},x^{2/3}\log^{2}x/h^{5/3})|+|S(x^{2/3}\log^{2}x/h^{5/3},2\sqrt{x})|.$

Thus, Theorem 2 is established by the following three propositions as they
would imply $|S|\leq\frac{11h}{12}$ and, hence, there must be some squarefree
number in $(x,x+h]$.

###### Proposition 1

For $x^{1/6}<h<x^{1/5}$ and some sufficiently small $c>0$, we have

$|S(h,ch^{5/2}/x^{1/4})|\leq\frac{h}{12}.$

###### Proposition 2

Let $c>0$ be a fixed small constant. Then, for $h=C_{0}x^{5/26}$ with some
sufficiently large $C_{0}>0$, we have

$|S(ch^{5/2}/x^{1/4},x^{2/3}\log^{2}x/h^{5/3})|\leq\frac{h}{12}.$

###### Proposition 3

For $h<x^{1/5}$, we have

$|S(x^{2/3}\log^{2}x/h^{5/3},2\sqrt{x})|\leq\frac{h}{12}.$

##  3 Strategies in establishing Propositions 1, 2 and 3

Let $\lambda=1+1/L$ for some big positive parameter $L$. The general strategy
in proving Propositions 1, 2 and 3 is to divide the intervals of consideration
for $n$ into $\lambda$-adic intervals $(N,\lambda N]$. Following Filaseta and
Trifonov [4], we consider consecutive elements in $S(N,\lambda N)$ and define

$T_{N,\lambda}(b):=\\{n:n\text{ and }n+b\text{ are consecutive elements in
}S(N,\lambda N)\\},\text{ and }t_{N,\lambda}(b):=|T_{N,\lambda}(b)|.$ (3)

First, note that

$|S(N,\lambda N)|\leq 1+\sum_{b_{1}=1}^{\infty}t_{N,\lambda}(b_{1}).$

Clearly,

$\frac{N}{L}\geq\sum_{b_{1}=1}^{\infty}b_{1}t_{N,\lambda}(b_{1})\geq\sum_{b_{1}>B}b_{1}t_{N,\lambda}(b_{1})\geq
B\sum_{b_{1}>B}t_{N,\lambda}(b_{1})$

which implies

$\sum_{b_{1}>B}t_{N,\lambda}(b_{1})\leq\frac{N}{LB}\leq\frac{h}{L\log^{1.2}x}$

for $B=N\log^{1.2}x/h$. By Lemma 3 below, we know that of every three
consecutive elements in $S(N,\lambda N)$, there are two consecutive elements
with distance $>R=cN^{4/3}/x^{1/3}$ from one another. Hence,

$\sum_{b_{1}\leq R}t_{N,\lambda}(b_{1})\leq
1+\sum_{b_{1}>R}t_{N,\lambda}(b_{1}),$

and

$|S(N,\lambda N)|\leq 2+2\sum_{b_{1}>R}t_{N,\lambda}(b_{1})\ll
1+\frac{h}{L\log^{1.2}x}+\sum_{R<b_{1}\leq
N\log^{1.2}x/h}t_{N,\lambda}(b_{1}).$ (4)

So, from now on, we restrict our attention to $R<b_{1}\leq N\log^{1.2}x/h$. In
[4], Filaseta and Trifonov established the bound

$t_{N,2}(b)\ll\frac{b^{1/3}x^{(3\theta+1)/3}\log
x}{N^{5}}+\frac{b^{3}x}{N^{4}}$ (5)

by considering certain differencings on the function $f(n)=x/n^{2}$. In a
subsequent paper [5], Filastea and Trifonov distillated their differencing
method and came up with the following general theorem on small fractional
parts:

###### Theorem 3

Let $N>1$. Let $r\geq 3$ be an integer. Let $T$ be a positive real number.
Suppose that $f$ is a function with at least $r$ derivatives with
$f^{(j)}(u)\asymp TN^{-j}$ for $j\in\\{r-2,r-1,r\\}$ and $u\in(N,2N]$. Let
$\delta$ be a positive real number with

$\delta<k\min\\{TN^{-r+2},T^{(r-4)/(r-2)}N^{-r+3}+TN^{-r+1}\\},$

for some sufficiently small constant $k$ depending on $r$ and the implied
constants in the asympotic formulas $f^{(j)}(u)\asymp TN^{-j}$ above. Let

$S=\\{u\in(N,2N]\cap\mathbb{Z}:||f(u)||<\delta\\}.$

Then

$|S|\ll
T^{\frac{2}{r(r+1)}}N^{\frac{r-1}{r+1}}+N\delta^{\frac{2}{(r-1)(r-2)}}+N(\delta
TN^{1-r})^{\frac{1}{r^{2}-3r+4}}.$

We will use this to establish Proposition 1. To prove Proposition 2, we will
apply the following proposition which allows one to extend the range of
validity for $S(N,\lambda N)=o(h)$ to be true.

###### Proposition 4

Suppose $N\leq C_{\alpha}x^{\alpha}$ for some $\theta<\alpha<2/5$ and some
constant $C_{\alpha}>0$. Then

$\sum_{Q<b_{1}\leq
N\log^{1.2}x/h}t_{N,\lambda}(b_{1})\ll\frac{N}{h}+\frac{x^{1/4}}{LN^{1/4}}+\frac{hx^{1/3}}{N^{7/5+4\delta/3}}+\frac{x}{N^{16/5+4\delta}}+\frac{hN^{4-2\beta}}{Lx^{1-\theta}}$
(6)

for

$C_{\alpha}^{\prime}\min(x^{\frac{1}{4-\beta}},x^{\frac{5}{23-9\beta}},x^{\frac{1}{27/5-3\beta-2\delta}})<N\leq
C_{\alpha}x^{\alpha},\;\;2-\frac{1-\theta}{2\alpha}<\beta<1$ (7)

and some sufficiently large constant $C_{\alpha}^{\prime}>0$. The lower bound
for $\beta$ in (7) is to guarantee that the last error term in (6) is much
less than $h$.

To prove Proposition 3, we are going to establish asymptotic relations between
shifts $a$ and $b$ when both $mn^{2}$ and $(m-a)(n+b)^{2}$ are members of
$S(N,\lambda N)$. Then we apply them to study quadruples
$(m+a_{1})(n-b_{1})^{2}\approx
mn^{2}\approx(m-a_{2})(n+b_{2})^{2}\approx(m-a_{2}-a_{3})(n+b_{2}+b_{3})^{2}$
and generalize Roth differencing and Filaseta-Trifonov differencing by
allowing $b_{1}$ to be different from $b_{3}$. Another new ingredient is the
introduction of a new differencing and using of the interplay among these
three differencings.

Here, let us give some comments on the propositions. Proposition 1 provides a
good estimate for small values of $N$ by extending it slightly beyond $N=h$
(Filaseta and Trifonov stopped at $S_{1}$ in (2) for the initial segment). On
the other end, Proposition 3 provides a $o(h)$ bound for $S(N,\lambda N)$ for
$N$ well below $x^{2/5}$. This is a substantial improvement, as previously,
this is only known to be true for $N>x^{2/5}$ when $h=x^{1/5}$. Filaseta and
Trifonov hinted in [4] that future improvement might depend on extending this
range. Finally, Proposition 2 provides a bridge for the middle range through
the use of the new differencing and a recursive argument via Proposition 4.

This paper is organized as follows. First, we will prove Proposition 1 via
Theorem 3. Then, we will study the asymptotic relations between shifts $a$ and
$b$ when $mn^{2}$ and $(m-a)(n+b)^{2}$ are members of $S(N,\lambda N)$.
Afterwards, we will revisit Roth differencing and Filaseta-Trifonov
differencing as introduced in [4] through these new asymptotic relations. We
will generalize these two differencings and introduce a third one. Then, we
will apply these to prove Proposition 3. Finally, we will prove Proposition 4
by considering consecutive elements $n,n+b_{2}\in T_{N,\lambda}(b_{1})$
depending on the size of the product $b_{1}b_{2}$. This will then be used to
establish Proposition 2.

##  4 Proof of Proposition 1

Recall $M_{n}$ denotes the number of multiples of $n^{2}$ in the interval
$(x,x+h]$. In this section, we will focus on

$|S(h,ch^{5/2}/x^{1/4})|=\sum_{h<n\leq ch^{5/2}/x^{1/4}}M_{n}=\sum_{h<n\leq
ch^{5/2}/x^{1/4}}\Bigl{[}\frac{x+h}{n^{2}}\Bigr{]}-\Bigl{[}\frac{x}{n^{2}}\Bigr{]}$
(8)

We divide $n$ into dyadic intervals, say $n\in(N,2N]$ with $h\leq N\leq
ch^{5/2}/x^{1/4}$. We will focus on bounding

$|S(N,2N)|=\sum_{N<n\leq 2N}M_{n}=\sum_{N<n\leq
2N}\Bigl{[}\frac{x+h}{n^{2}}\Bigr{]}-\Bigl{[}\frac{x}{n^{2}}\Bigr{]}.$ (9)

By the discussion in section 2, the summand in (9) is $1$ if and only if there
is some integer $m$ such that

$\frac{x}{n^{2}}<m\leq\frac{x+h}{n^{2}},$ (10)

and $0$ otherwise. One can see that if condition (10) holds for some
$n\in(N,2N]$, then $\|\frac{x}{n^{2}}\|<\frac{h}{N^{2}}$. Thus,

$|S(N,2N)|\leq\Big{|}\bigl{\\{}n\in(N,2N]:\big{\|}\frac{x}{n^{2}}\big{\|}<\frac{h}{N^{2}}\big{\\}}\Big{|}=:M_{N}.$
(11)

Let

$f(u)=\frac{x}{u^{2}},\;\;\delta=\frac{h}{N^{2}},\;\;\text{ and
}\;\;T=\frac{x}{N^{2}}.$

Now, we apply Theorem 3 with $r=4$. One can check that the conditions in
Theorem 3 are satisfied for $h\leq N\leq x^{1/4}$ and $h<x^{1/5}$. Therefore,

$M_{N}\ll T^{1/10}N^{3/5}+N\delta^{1/3}+N\Bigl{(}\frac{\delta
T}{N^{3}}\Bigr{)}^{1/8}\ll x^{1/10}N^{2/5}$ (12)

as the first term dominates due to $x^{1/6}<h<x^{1/5}$. Summing (12) over the
dyadic intervals, we obtain

$|S(h,ch^{5/2}/x^{1/4})|<h/12.$

for sufficiently small $c>0$. This establishes Proposition 1.

##  5 Asymptotic relations between shifts

Suppose $mn^{2}$ and $(m-a)(n+b)^{2}$ are members in $S$ for some integers
$1\leq a<m/2$, $1\leq b<n/2$ and $n\geq h$. Consider

$d:=(m-a)(n+b)^{2}-mn^{2}.$ (13)

Then $d=(2bm-an)n+mb^{2}-2abn-ab^{2}$ which gives

$\Delta:=2bm-
an=2ab-\frac{mb^{2}}{n}+\frac{ab^{2}}{n}+\frac{d}{n}=\frac{3}{2}ab+\frac{(an-2mb)b}{2n}+\frac{ab^{2}}{n}+\frac{d}{n}.$

Hence,

$\Bigl{(}1+\frac{b}{2n}\Bigr{)}\Delta=\frac{3}{2}ab+\frac{ab^{2}}{n}+\frac{d}{n}$

which gives

$\Delta=\frac{3ab}{2(1+\frac{b}{2n})}+\frac{ab^{2}}{n(1+\frac{b}{2n})}+\frac{d}{n(1+\frac{b}{2n})}$
(14)

In particular, $ab\leq\Delta\leq 3ab$ for $x$ (and hence $n$) sufficiently
large. From (13), we also have $\frac{d}{(n+b)^{2}}=m-a-m(\frac{n}{n+b})^{2}$
which gives

$a=m\Bigl{[}1-\frac{1}{(1+\frac{b}{n})^{2}}\Bigr{]}-\frac{d}{(n+b)^{2}}=\frac{2mb}{n}\cdot\frac{1+\frac{b}{2n}}{(1+\frac{b}{n})^{2}}-\frac{d}{(n+b)^{2}}.$
(15)

Using the geometric series formula
$\frac{1-x^{k}}{1-x}=1+x+x^{2}+...+x^{k-1}$, one obtains

$a=\frac{2mb}{n}-\frac{3mb^{2}}{n^{2}}+\frac{4mb^{3}}{n^{3}}-\frac{5mb^{4}}{n^{4}}+\frac{6mb^{5}}{n^{5}}+O\Bigl{(}\frac{mb^{6}}{n^{6}}\Bigr{)}-\frac{d}{(n+b)^{2}}.$
(16)

Summarizing all these, we have the following lemma.

###### Lemma 1

Suppose $mn^{2}$ and $(m-a)(n+b)^{2}$ are members in $S$ for some integers
$n\geq h$ and $1\leq b<n/2$. Then $a$ and $b$ satisfy the following relations:

$a=\frac{2mb}{n}-\frac{3mb^{2}}{n^{2}}+\frac{4mb^{3}}{n^{3}}-\frac{5mb^{4}}{n^{4}}+\frac{6mb^{5}}{n^{5}}+O\Bigl{(}\frac{mb^{6}}{n^{6}}\Bigr{)}-\frac{d}{(n+b)^{2}},$

and

$\Delta:=2bm-
an=\frac{3ab}{2(1+\frac{b}{2n})}+\frac{ab^{2}}{n(1+\frac{b}{2n})}+\frac{d}{n(1+\frac{b}{2n})}$

where $d=(m-a)(n+b)^{2}-mn^{2}$.

A similar calculation yields

###### Lemma 2

Suppse $mn^{2}$ and $(m+a)(n-b)^{2}$ are members in $S$ for some integers
$n\geq h$ and $1\leq b<n/2$. Then $a$ and $b$ satisfy the following relations:

$a=\frac{2mb}{n}+\frac{3mb^{2}}{n^{2}}+\frac{4mb^{3}}{n^{3}}+\frac{5mb^{4}}{n^{4}}+O\Bigl{(}\frac{mb^{5}}{n^{5}}\Bigr{)}+\frac{d^{\prime}}{(n-b)^{2}},$

and

$\Delta:=2bm-
an=-\frac{3ab}{2(1-\frac{b}{2n})}+\frac{ab^{2}}{n(1-\frac{b}{2n})}-\frac{d^{\prime}}{n(1-\frac{b}{2n})}$

where $d^{\prime}=(m+a)(n-b)^{2}-mn^{2}$.

The idea of using differencing was introduced by Roth [10] and is the building
block for subsequent improvements via elementary method. We are going to
illustrate how one can combine differencing and our asymptotic relations in
Lemmas 1 and 2 to recover the following result used in previous works (cf.
[3]):

###### Lemma 3

If $I$ is a subinterval of $(N,\lambda N]$ with $N\geq h$ and
$|I|\leq\frac{c_{0}N^{4/3}}{x^{1/3}}=:R$ with some small $0<c_{0}<1$, then
$|S(N,\lambda N)\cap I|\leq 2$.

Proof: Suppose $(m+a_{1})(n-b_{1})^{2}$, $mn^{2}$ and $(m-a_{2})(n+b_{2})^{2}$
are members in $S(N,\lambda N)$ for some integers $1\leq b_{1},b_{2}<N/L$ and
$h\leq n$. We consider the difference

$D:=b_{1}(m-a_{2})(n+b_{2})^{2}-(b_{1}+b_{2})mn^{2}+b_{2}(m+a_{1})(n-b_{1})^{2}.$

Also, denote

$\Delta_{1}=2b_{1}m-a_{1}n,\;\text{ and }\;\Delta_{2}=2b_{2}m-a_{2}n.$

By some simple algebra, one gets

$D=b_{1}b_{2}(b_{1}+b_{2})m+(a_{1}b_{2}-a_{2}b_{1})n^{2}-2b_{1}b_{2}(a_{1}+a_{2})n+b_{1}b_{2}(a_{1}b_{1}-a_{2}b_{2}).$

Substituting $a_{1}=\frac{2b_{1}m}{n}-\frac{\Delta_{1}}{n}$ and
$a_{2}=\frac{2b_{2}m}{n}-\frac{\Delta_{2}}{n}$ into the third term above, we
have

$D=(a_{1}b_{2}-a_{2}b_{1})n^{2}-3b_{1}b_{2}(b_{1}+b_{2})m+2b_{1}b_{2}(\Delta_{1}+\Delta_{2})+b_{1}b_{2}(a_{1}b_{1}-a_{2}b_{2}).$
(17)

By the formulas of $\Delta$ in Lemmas 1 and 2,

$2b_{1}b_{2}(\Delta_{1}+\Delta_{2})+b_{1}b_{2}(a_{1}b_{1}-a_{2}b_{2})=-2b_{1}b_{2}(a_{1}b_{1}-a_{2}b_{2})+\frac{b_{1}b_{2}(a_{1}b_{1}^{2}+a_{2}b_{2}^{2})}{2n}+O\Bigl{(}\frac{b_{1}b_{2}(a_{1}b_{1}^{3}+a_{2}b_{2}^{3})}{n^{2}}+\frac{b_{1}b_{2}h}{n}\Bigr{)}$
(18)

which has size much smaller than the term $3b_{1}b_{2}(b_{1}+b_{2})m$ by the
formulas of $a$ in Lemmas 1 and 2 as $b_{1},b_{2}<N/L$. Putting (18) into
(17), one has

$D=(a_{1}b_{2}-a_{2}b_{1})n^{2}-3b_{1}b_{2}(b_{1}+b_{2})m-2b_{1}b_{2}(a_{1}b_{1}-a_{2}b_{2})+\frac{b_{1}b_{2}(a_{1}b_{1}^{2}+a_{2}b_{2}^{2})}{2n}+O\Bigl{(}\frac{b_{1}b_{2}(a_{1}b_{1}^{3}+a_{2}b_{2}^{3})}{n^{2}}+\frac{b_{1}b_{2}h}{n}\Bigr{)}.$
(19)

As $|D|\leq(b_{1}+b_{2})h\leq 2n^{2}/L$,

$a_{1}b_{2}-a_{2}b_{1}=\frac{3mb_{1}b_{2}(b_{1}+b_{2})}{n^{2}}+\frac{2b_{1}b_{2}(a_{1}b_{1}-a_{2}b_{2})}{n^{2}}-\frac{b_{1}b_{2}(a_{1}b_{1}^{2}+a_{2}b_{2}^{2})}{2n^{3}}+O\Bigl{(}\frac{b_{1}b_{2}(a_{1}b_{1}^{3}+a_{2}b_{2}^{3})}{n^{4}}+\frac{(b_{1}+b_{2})h}{n^{2}}\Big{)}$
(20)

If $a_{1}b_{2}-a_{2}b_{1}=0$, then (20) implies

$\frac{mb_{1}b_{2}(b_{1}+b_{2})}{n^{2}}\ll\frac{(b_{1}+b_{2})h}{n^{2}}\;\;\text{
or }\;\;mb_{1}b_{2}\ll h$

as $b_{1},b_{2}<N/L$ and $a_{1},a_{2}<10M/L$. If $m\geq x^{1/3}$, then
$x^{1/3}\leq mb_{1}b_{2}\ll h\ll x^{1/5}$ which is a contradiction. If $m\leq
x^{1/3}$, then $n\gg\sqrt{x/m}\geq x^{1/3}\geq m$. However, since
$a_{1},a_{2}$ are positive integers, Lemmas 1 and 2 imply that $b_{1},b_{2}\gg
n/m$. Hence $x^{1/3}\ll n\ll n^{2}/m=m(n/m)^{2}\ll mb_{1}b_{2}\ll h\ll
x^{1/5}$ which is impossible. Therefore, $a_{2}b_{1}-a_{1}b_{2}\neq 0$ and
(20) can be rewritten as

$a_{1}b_{2}-a_{2}b_{1}+O\Bigl{(}\frac{(b_{1}+b_{2})h}{n^{2}}\Big{)}=\frac{3mb_{1}b_{2}(b_{1}+b_{2})}{n^{2}}+O\Bigl{(}\frac{b_{1}b_{2}(a_{1}b_{1}+a_{2}b_{2})}{n^{2}}\Big{)}.$

One can easily see that the error term on each side is much smaller than the
main term on each side. Hence, $a_{1}b_{2}-a_{2}b_{1}\geq 1$ as the right hand
side is clearly positive. Thus,

$n^{2}\leq(a_{1}b_{2}-a_{2}b_{1})n^{2}\ll b_{1}b_{2}(b_{1}+b_{2})m$ (21)

which implies

$\max(b_{1},b_{2})\gg\Bigl{(}\frac{n^{2}}{m}\Bigr{)}^{1/3}\gg\frac{N^{4/3}}{x^{1/3}}$
(22)

and, hence, the lemma by picking $c_{0}$ sufficiently small.

##  6 Roth, Filaseta-Trifonov and a new differencing

Suppose $mn^{2},(m-a)(n+b)^{2}\in S(N,\lambda N)$, Roth [10] considered the
difference

$D_{1}:=m(2n-b)-(m-a)(2n+3b)=2an-4bm+3ab.$

By Lemma 1,

$\displaystyle D_{1}=$ $\displaystyle
2n\Bigl{[}\frac{2mb}{n}-\frac{3mb^{2}}{n^{2}}+\frac{4mb^{3}}{n^{3}}-\frac{5mb^{4}}{n^{4}}+\frac{6mb^{5}}{n^{5}}\Bigr{]}-4bm+3b\Bigl{[}\frac{2mb}{n}-\frac{3mb^{2}}{n^{2}}+\frac{4mb^{3}}{n^{3}}-\frac{5mb^{4}}{n^{4}}\Bigr{]}$
$\displaystyle+O\Bigl{(}\frac{h}{n}+\frac{mb^{6}}{n^{5}}\Bigr{)}$
$\displaystyle=$
$\displaystyle-\frac{mb^{3}}{n^{2}}+\frac{2mb^{4}}{n^{3}}-\frac{3mb^{5}}{n^{4}}+O\Bigl{(}\frac{h}{n}+\frac{mb^{6}}{n^{5}}\Bigr{)}$
(23)

as initial terms cancel out. This and the quantity $\Delta$ in Lemma 1 might
indicate how Roth came up with such a differencing. Similarly, suppose
$(m+a)(n-b)^{2},mn^{2}\in S(N,\lambda N)$. Then

$\displaystyle D_{1}^{\prime}=$
$\displaystyle(m+a)(2n-3b)-m(2n+b)=2an-4bm-3ab$ $\displaystyle=$
$\displaystyle-\frac{mb^{3}}{n^{2}}-\frac{2mb^{4}}{n^{3}}-\frac{3mb^{5}}{n^{4}}+O\Bigl{(}\frac{h}{n}+\frac{mb^{6}}{n^{5}}\Bigr{)}$
(24)

Now, suppose $(m+a_{1})(n-b_{1})^{2}$, $mn^{2}$, $(m-a_{2})(n+b_{2})^{2}$ and
$(m-a_{2}-a_{3})(n+b_{2}+b_{3})^{2}$ are four members in $S(N,\lambda N)$ for
some integers $R<b_{1},b_{2},b_{3}$ and $h\leq n$. Clearly,
$b_{1},b_{2},b_{3}<N/L$ as $\lambda=1+1/L$. Applying (6) to
$(m+a_{1})(n-b_{1})^{2}$, $mn^{2}$ and (6) to $(m-a_{2})(n+b_{2})^{2}$,
$(m-a_{2}-a_{3})(n+b_{2}+b_{3})^{2}$, we obtain

$\displaystyle D_{2}:=$
$\displaystyle(m+a_{1})(2n-3b_{1})-m(2n+b_{1})=2a_{1}n-4b_{1}m-3a_{1}b_{1}$
$\displaystyle=$
$\displaystyle-\frac{mb_{1}^{3}}{n^{2}}-\frac{2mb_{1}^{4}}{n^{3}}-\frac{3mb_{1}^{5}}{n^{4}}+O\Bigl{(}\frac{h}{n}+\frac{mb_{1}^{6}}{n^{5}}\Bigr{)}$
(25)

and

$\displaystyle D_{3}:=$
$\displaystyle(m-a_{2})(2(n+b_{2})-b_{3})-(m-a_{2}-a_{3})(2(n+b_{2})+3b_{3})=2a_{3}n-4b_{3}m+2a_{3}b_{2}+4a_{2}b_{3}+3a_{3}b_{3}$
$\displaystyle=$
$\displaystyle-\frac{(m-a_{2})b_{3}^{3}}{(n+b_{2})^{2}}+\frac{2(m-a_{2})b_{3}^{4}}{(n+b_{2})^{3}}-\frac{3(m-a_{2})b_{3}^{5}}{(n+b_{2})^{4}}+O\Bigl{(}\frac{h}{n}+\frac{mb_{3}^{6}}{n^{5}}\Bigr{)}.$
(26)

Subtracting (6) from (6), we have

$\displaystyle D_{4}:=$ $\displaystyle D_{3}-D_{2}$ $\displaystyle=$
$\displaystyle
2(a_{3}-a_{1})n-4(b_{3}-b_{1})m+2a_{3}b_{2}+4a_{2}b_{3}+3a_{1}b_{1}+3a_{3}b_{3}$
$\displaystyle=$
$\displaystyle\frac{mb_{1}^{3}}{n^{2}}-\frac{mb_{3}^{3}}{n^{2}}+\Bigl{[}\frac{mb_{3}^{3}}{n^{2}}-\frac{(m-a_{2})b_{3}^{3}}{n^{2}}\Bigr{]}+\Bigl{[}\frac{(m-a_{2})b_{3}^{3}}{n^{2}}-\frac{(m-a_{2})b_{3}^{3}}{(n+b_{2})^{2}}\Bigr{]}$
$\displaystyle+\frac{2mb_{1}^{4}}{n^{3}}+\frac{2mb_{3}^{4}}{n^{3}}-\Bigl{[}\frac{2mb_{3}^{4}}{n^{3}}-\frac{2(m-a_{2})b_{3}^{4}}{n^{3}}\Bigr{]}-\Bigl{[}\frac{2(m-a_{2})b_{3}^{4}}{n^{3}}-\frac{2(m-a_{2})b_{3}^{4}}{(n+b_{2})^{3}}\Bigr{]}$
$\displaystyle+\frac{3mb_{1}^{5}}{n^{4}}-\frac{3mb_{3}^{5}}{n^{4}}+\Bigl{[}\frac{3mb_{3}^{5}}{n^{4}}-\frac{3(m-a_{2})b_{3}^{5}}{n^{4}}\Bigr{]}+\Bigl{[}\frac{3(m-a_{2})b_{3}^{5}}{n^{4}}-\frac{3(m-a_{2})b_{3}^{5}}{(n+b_{2})^{4}}\Bigr{]}+O\Bigl{(}\frac{h}{n}+\frac{m(b_{1}^{6}+b_{3}^{6})}{n^{5}}\Bigr{)}$
$\displaystyle=$
$\displaystyle\frac{m(2b_{1}^{4}+4b_{2}b_{3}^{3}+2b_{3}^{4})}{n^{3}}+\frac{m(b_{1}-b_{3})(b_{1}^{2}+b_{1}b_{3}+b_{3}^{2})}{n^{2}}$
$\displaystyle+\frac{m(3b_{1}^{5}-10b_{2}^{2}b_{3}^{3}-10b_{2}b_{3}^{4}-3b_{3}^{5})}{n^{4}}+O\Bigl{(}\frac{h}{n}+\frac{m(b_{1}^{3}+b_{3}^{3})(b_{1}+b_{2}+b_{3})^{3}}{n^{5}}\Bigr{)}$
(27)

by Lemma 1. This differencing, $D_{4}$, was introduced and used in [4] and we
call it Roth differencing.

To supplement Roth differencing, Filaseta and Trifonov [4, section 3] also
introduced the following second differencing (with slightly different
notations):

$\displaystyle D^{\prime}:=$ $\displaystyle
b_{2}(m+a_{1})-(b_{1}+b_{2}+b_{3})m+(b_{1}+b_{2}+b_{3})(m-a_{2})-b_{2}(m-a_{2}-a_{3})$
$\displaystyle=$
$\displaystyle(b_{2}a_{1}-b_{1}a_{2})-(b_{3}a_{2}-b_{2}a_{3}).$

It is more general here as we do not require $b_{1}=b_{3}$. In any case, we
call this Filaseta-Trifonov differencing. By (20) in Lemma 3,

$\displaystyle D^{\prime}=$
$\displaystyle\frac{3mb_{1}b_{2}(b_{1}+b_{2})}{n^{2}}-\frac{3(m-a_{2})b_{2}b_{3}(b_{2}+b_{3})}{(n+b_{2})^{2}}+\frac{2b_{1}b_{2}(a_{1}b_{1}-a_{2}b_{2})}{n^{2}}-\frac{2b_{2}b_{3}(a_{2}b_{2}-a_{3}b_{3})}{(n+b_{2})^{2}}$
$\displaystyle-\frac{b_{1}b_{2}(a_{1}b_{1}^{2}+a_{2}b_{2}^{2})}{2n^{3}}+\frac{b_{2}b_{3}(a_{2}b_{2}^{2}+a_{3}b_{3}^{2})}{2(n+b_{2})^{3}}+O\Bigl{(}\frac{b_{2}(b_{1}+b_{3})(a_{1}b_{1}^{3}+a_{2}b_{2}^{3}+a_{3}b_{3}^{3})}{n^{4}}+\frac{(b_{1}+b_{2}+b_{3})h}{n^{2}}\Bigr{)}.$
(28)

By Lemmas 1 and 2, the first two terms in (6) can be rewritten as

$\displaystyle\Bigl{[}\frac{3mb_{1}b_{2}(b_{1}+b_{2})}{n^{2}}-\frac{3mb_{3}b_{2}(b_{1}+b_{2})}{n^{2}}\Bigr{]}+\Bigl{[}\frac{3mb_{3}b_{2}(b_{1}+b_{2})}{n^{2}}-\frac{3mb_{3}b_{2}(b_{3}+b_{2})}{n^{2}}\Bigr{]}$
$\displaystyle+\Bigl{[}\frac{3mb_{2}b_{3}(b_{2}+b_{3})}{n^{2}}-\frac{3(m-a_{2})b_{2}b_{3}(b_{2}+b_{3})}{n^{2}}\Bigr{]}+\Bigl{[}\frac{3(m-a_{2})b_{2}b_{3}(b_{2}+b_{3})}{n^{2}}-\frac{3(m-a_{2})b_{2}b_{3}(b_{2}+b_{3})}{(n+b_{2})^{2}}\Bigr{]}$
$\displaystyle=$
$\displaystyle\frac{3m(b_{1}-b_{3})b_{2}(b_{1}+b_{2}+b_{3})}{n^{2}}+\frac{12mb_{2}^{2}b_{3}(b_{2}+b_{3})}{n^{3}}-\frac{30mb_{2}^{3}b_{3}(b_{2}+b_{3})}{n^{4}}+O\Bigl{(}\frac{mb_{2}^{4}b_{3}(b_{2}+b_{3})}{n^{5}}+\frac{(b_{2}+b_{3})h}{n^{2}}\Bigr{)},$
(29)

the third and fourth terms in (6) can be rewritten as

$\displaystyle\frac{2a_{1}b_{1}^{2}b_{2}-2a_{2}b_{2}^{2}(b_{1}+b_{3})+2a_{3}b_{3}^{2}b_{2}}{n^{2}}+\frac{4b_{2}^{2}b_{3}(a_{2}b_{2}-a_{3}b_{3})}{n^{3}}+O\Bigl{(}\frac{b_{2}^{3}b_{3}(a_{2}b_{2}-a_{3}b_{3})}{n^{4}}\Bigr{)}$
$\displaystyle=$
$\displaystyle\frac{4m(b_{1}^{3}b_{2}-b_{2}^{3}(b_{1}+b_{3})+b_{3}^{3}b_{2})}{n^{3}}+\frac{m(6b_{1}^{4}b_{2}+6b_{2}^{4}b_{1}+14b_{2}^{4}b_{3}-28b_{2}^{2}b_{3}^{3}-6b_{2}b_{3}^{4})}{n^{4}}$
(30)
$\displaystyle+O\Bigl{(}\frac{m(b_{1}+b_{3})(b_{1}^{5}+b_{2}^{5}+b_{3}^{5})}{n^{5}}+\frac{(b_{1}+b_{2}+b_{3})h}{n^{2}}\Bigr{)},$

and the fifth and sixth terms in (6) can be rewritten as

$\displaystyle\frac{-a_{1}b_{1}^{3}b_{2}-a_{2}b_{1}b_{2}^{3}+a_{2}b_{3}b_{2}^{3}+a_{3}b_{3}^{3}b_{2}}{2n^{3}}+O\Bigl{(}\frac{b_{2}^{2}(b_{1}+b_{3})(a_{1}b_{1}^{2}+a_{2}b_{2}^{2}+a_{3}b_{3}^{2})}{n^{4}}\Bigr{)}$
$\displaystyle=$
$\displaystyle\frac{-b_{1}^{4}b_{2}-b_{1}b_{2}^{4}+b_{3}b_{2}^{4}+b_{3}^{4}b_{2}}{n^{4}}+O\Bigl{(}\frac{mb_{2}^{2}(b_{1}+b_{3})(b_{1}^{3}+b_{2}^{3}+b_{3}^{3})}{n^{5}}\Bigr{)}$
(31)

Putting (6), (6) and (6) into (6), we get

$\displaystyle D^{\prime}=$
$\displaystyle\frac{4mb_{2}(b_{1}^{3}+b_{3}^{3}+3b_{3}^{2}b_{2}+2b_{3}b_{2}^{2}-b_{1}b_{2}^{2})}{n^{3}}+\frac{3m(b_{1}-b_{3})b_{2}(b_{1}+b_{2}+b_{3})}{n^{2}}$
$\displaystyle+\frac{m(5b_{1}^{4}b_{2}+5b_{2}^{4}b_{1}-15b_{2}^{4}b_{3}-30b_{2}^{3}b_{3}^{2}-28b_{2}^{2}b_{3}^{3}-5b_{2}b_{3}^{4})}{n^{4}}$
(32)
$\displaystyle+O\Bigl{(}\frac{m(b_{1}+b_{3})(b_{1}^{5}+b_{2}^{5}+b_{3}^{5})}{n^{5}}+\frac{(b_{1}+b_{2}+b_{3})h}{n^{2}}\Bigr{)}.$

For $h=Cx^{5/26}$ and $n\gg h^{5/2}/x^{1/4}\gg h^{6/5}$, it is worth
mentioning that when $b_{1}=b_{3}$ and $b_{2}\gg b_{1}$,

$D^{\prime}=\frac{4mb_{1}b_{2}(b_{2}+b_{1})(b_{2}+2b_{1})}{n^{3}}+O\Bigl{(}\frac{mb_{1}(b_{1}^{4}+b_{2}^{4})}{n^{4}}+\frac{(b_{1}+b_{2})h}{n^{2}}\Bigr{)}>0$
(33)

as the main term is much bigger than the error terms by $b_{1}+b_{2}<n/L$ and
(21).

Now, we introduce a new differencing, namely:

$a_{1}-a_{3}.$

Note that $b_{1}+b_{2}+b_{3}<N/L$. By Lemmas 1 and 2,

$a_{1}=\frac{2mb_{1}}{n}+\frac{3mb_{1}^{2}}{n^{2}}+O\Bigl{(}\frac{h}{n^{2}}+\frac{mb_{1}^{3}}{n^{3}}\Bigr{)}$
(34)

and

$a_{3}=\frac{2(m-a_{2})b_{3}}{n+b_{2}}-\frac{3(m-a_{2})b_{3}^{2}}{(n+b_{2})^{2}}+O\Bigl{(}\frac{h}{n^{2}}+\frac{mb_{3}^{3}}{n^{3}}\Bigr{)}.$
(35)

As $b_{1},b_{2},b_{3}<N/L$, Lemma 1 implies $0<a_{1},a_{2},a_{3}<10M/L$.
Subtracting (35) from (34), we obtain

$\displaystyle a_{1}-a_{3}=$
$\displaystyle\Bigl{[}\frac{2mb_{1}}{n}-\frac{2(m-a_{2})b_{1}}{n}\Bigr{]}+\Bigl{[}\frac{2(m-a_{2})b_{1}}{n}-\frac{2(m-a_{2})b_{1}}{n+b_{2}}\Bigr{]}+\Bigl{[}\frac{2(m-a_{2})b_{1}}{n+b_{2}}-\frac{2(m-a_{2})b_{3}}{n+b_{2}}\Bigr{]}$
$\displaystyle+\frac{3mb_{1}^{2}}{n^{2}}+\frac{3(m-a_{2})b_{1}^{2}}{(n+b_{2})^{2}}-\Bigl{[}\frac{3(m-a_{2})b_{1}^{2}}{(n+b_{2})^{2}}-\frac{3(m-a_{2})b_{3}^{2}}{(n+b_{2})^{2}}\Bigr{]}+O\Bigl{(}\frac{h}{n^{2}}+\frac{m(b_{1}^{3}+b_{3}^{3})}{n^{3}}\Bigr{)}$
$\displaystyle=$
$\displaystyle\frac{6mb_{1}(b_{1}+b_{2})}{n^{2}}+\frac{2m(b_{1}-b_{3})}{n}+O\Bigl{(}\frac{h}{n^{2}}+\frac{m(b_{1}+b_{3})(b_{1}+b_{2}+b_{3})^{2}}{n^{3}}+\frac{m|b_{1}-b_{3}|(b_{1}+b_{2}+b_{3})}{n^{2}}\Bigr{)}$
(36)

by Lemma 1 and (20). Equation (34) also implies that $\frac{mb_{1}}{n}\gg
a_{1}\geq 1$ or $b_{1}\gg\frac{n}{m}$.

Case 1: $N\leq x^{2/5}$. We use the bound $b_{1}>R$ and get

$\frac{mb_{1}(b_{1}+b_{2})}{n^{2}}\gg\frac{m}{n^{2}}b_{1}^{2}\gg\frac{x}{N^{4}}\Bigl{(}\frac{N^{4/3}}{x^{1/3}}\Bigr{)}^{2}=\frac{x^{1/3}}{N^{4/3}}\geq\frac{x^{1/5}}{N}.$

Case 2: $N>x^{2/5}$. We use the bound $b_{1}\gg\frac{n}{m}\gg\frac{N^{3}}{x}$
and get

$\frac{mb_{1}(b_{1}+b_{2})}{n^{2}}\gg\frac{m}{n^{2}}b_{1}^{2}\gg\frac{N^{2}}{x}>\frac{x^{1/5}}{N}.$

In either cases, the main term $\frac{6mb_{1}(b_{1}+b_{2})}{n^{2}}$ is much
bigger than the error term $O(\frac{h}{n^{2}})$ while the term
$\frac{2m(b_{1}-b_{3})}{n}$ is much bigger than the error term
$O(\frac{m|b_{1}-b_{3}|(b_{1}+b_{2}+b_{3})}{n^{2}})$ unless $b_{1}=b_{3}$ (in
which case both terms vanish).

##  7 Proof of Proposition 3

Here we consider $x^{2/3}\log^{2}x/h^{5/3}<n\leq 2\sqrt{x}$. We break this
interval into intervals of the form $(N,\lambda N]$ with
$x^{2/3}\log^{2}x/h^{5/3}\leq N\leq 2\sqrt{x}$. Note that, as $h<x^{1/5}$, we
have $N\geq x^{1/3}\log^{2}x$ and $M\ll x^{1/3}/\log^{4}x$. Suppose
$(m+a_{1})(n-b_{1})^{2}$, $mn^{2}$, $(m-a_{2})(n+b_{2})^{2}$ and
$(m-a_{2}-a_{3})(n+b_{2}+b_{3})^{2}$ are four members in $S(N,\lambda N)$ with
$n-b_{1}\in T_{N,\lambda}(b_{1})$, $n+b_{2}\in T_{N,\lambda}(b_{3})$ for some
integers $R<b_{1},b_{3}$ and

$|b_{1}-b_{3}|\leq\frac{0.0001N}{M}.$ (37)

Recall (6)

$a_{1}-a_{3}=\frac{6m(b_{1}+b_{2})b_{1}}{n^{2}}+\frac{2m(b_{1}-b_{3})}{n}+O\Bigl{(}\frac{h}{n^{2}}+\frac{m(b_{1}+b_{3})(b_{1}+b_{2}+b_{3})^{2}}{n^{3}}+\frac{m|b_{1}-b_{3}|(b_{1}+b_{2}+b_{3})}{n^{2}}\Bigr{)}.$
(38)

If $a_{1}-a_{3}=0$, then we have

$b_{1}-b_{3}+O\Bigl{(}\frac{|b_{1}-b_{3}|(b_{1}+b_{2}+b_{3})}{n}\Bigr{)}=-\frac{3(b_{1}+b_{2})b_{1}}{n}+O\Bigl{(}\frac{h}{mn}+\frac{(b_{1}+b_{3})(b_{1}+b_{2}+b_{3})^{2}}{n^{2}}\Bigr{)}$
(39)

By (21), the error terms on the right hand side are much less than the
corresponding main term. Thus, $b_{1}-b_{3}$ must be negative and
$b_{3}=b_{1}(1+O(1/L))$ as $b_{1}+b_{2}+b_{3}<n/L$. Putting (39) into (6) and
(6), we have

$D^{\prime}=-\frac{5mb_{1}b_{2}(b_{2}+b_{1})(b_{2}+2b_{1})}{n^{3}}\Bigl{(}1+O\Bigl{(}\frac{1}{L}\Bigr{)}\Bigr{)}+O\Bigl{(}\frac{(b_{1}+b_{2})h}{n^{2}}\Bigr{)}$
(40)

and

$D_{4}=-\frac{5mb_{1}^{3}(b_{1}+b_{2})}{n^{3}}\Bigl{(}1+O\Bigl{(}\frac{1}{L}\Bigr{)}\Bigr{)}+O\Bigl{(}\frac{h}{n}\Bigr{)}.$
(41)

This implies $|D^{\prime}|\geq 1$ by (21), and $D_{4}$ has size
$mb_{1}^{3}(b_{1}+b_{2})/n^{3}$ plus an error of $O(h/n)$. If $a_{1}-a_{3}\neq
0$, then (38) and condition (37) imply

$(b_{1}+b_{2})b_{1}\geq\frac{n^{2}}{9m}\;\;\text{ and
}\;\;|b_{1}-b_{3}|\leq\frac{0.0001N}{M}\leq\frac{0.001b_{1}(b_{1}+b_{2})}{n}\ll\frac{b_{1}}{L}.$

Putting this into (6) and (6), we have

$D^{\prime}\geq\frac{3.999mb_{1}b_{2}(b_{2}+b_{1})(b_{2}+2b_{1})}{n^{3}}\Bigl{(}1+O\Bigl{(}\frac{1}{L}\Bigr{)}\Bigr{)}+O\Bigl{(}\frac{(b_{1}+b_{2})h}{n^{2}}\Bigr{)}$
(42)

and

$D_{4}\geq\frac{3.999mb_{1}^{3}(b_{1}+b_{2})}{n^{3}}\Bigl{(}1+O\Bigl{(}\frac{1}{L}\Bigr{)}\Bigr{)}+O\Bigl{(}\frac{h}{n}\Bigr{)}.$
(43)

Again, this implies $|D^{\prime}|\geq 1$ by (21), and $D_{4}$ has size
$mb_{1}^{3}(b_{1}+b_{2})/n^{3}$ plus an error of $O(h/n)$.

Suppose $b_{1}\in(B,2B]$. We break this dyadic interval $(B,2B]$ into
subintervals $I$ of length $0.0001N/M$ according to condition (37). We want to
bound $\sum_{b_{1}\in I}t_{N,\lambda}(b_{1})$. Following [4], we break
$(N,\lambda N]$ into subintervals $J$ of length $\frac{cN^{3}}{MB^{3}}$ and
consider the set $U=(\cup_{b_{1}\in I}T_{N,\lambda}(b_{1}))\cap J$. Define

$U_{N,\lambda,I,J}(b_{2}):=\\{n:n\text{ and }n+b_{2}\text{ are consecutive
elements in }U\\},\text{ and
}u_{N,\lambda,I,J}(b_{2}):=|U_{N,\lambda,I,J}(b_{2})|.$

Clearly, one has

$\sum_{b_{1}\in
I}t_{N,\lambda}(b_{1})\leq\sum_{J}\Bigl{(}1+\sum_{b_{2}}u_{N,\lambda,I,J}(b_{2})\Bigr{)}.$
(44)

Say, for some $b_{1}\in I$, $n-b_{1}\in T_{N,\lambda}(b_{1})\cap J$ is an
element in $U_{N,\lambda,I,J}(b_{2})$ for some $b_{2}$. Then $n+b_{2}\in
T_{N,\lambda}(b_{3})\cap J$ for some $b_{3}\in I$. We must have $b_{2}>0$,
$b_{1}+b_{2}\leq\frac{cN^{3}}{MB^{3}}$ and $D_{4}=0$ if $c>0$ is sufficiently
small. The contribution from those $b_{2}\leq C^{\prime}B$ in (44) is small
because (40) and (42) imply

$1\leq|D^{\prime}|\ll\frac{Mb_{1}^{4}}{N^{3}}\;\;\text{ or
}\;\;b_{1}\gg\frac{N^{3/4}}{M^{1/4}},$

and

$\sum_{B}\sum_{I}\sum_{J}\sum_{b_{2}\leq
C^{\prime}B}u_{N,\lambda,I,J}(b_{2})\leq\sum_{b_{1}\gg
N^{3/4}/M^{1/4}}t_{N,\lambda}(b_{1})\ll\frac{N/L}{N^{3/4}/M^{1/4}}=\frac{x^{1/4}}{LN^{1/4}}.$
(45)

For $b_{2}>C^{\prime}b_{1}$ with large enough $C^{\prime}>0$, $D_{4}=0$, (41)
and (43) give $b_{1}+b_{2}\ll\frac{hN^{2}}{MB^{3}}$ while (40) and (42) imply
$b_{2}\gg b_{1}+b_{2}\gg\frac{N}{M^{1/3}B^{1/3}}$. These give upper and lower
bounds for the distance $b_{1}+b_{2}$ on any such quadruples and is a method
introduced in [4]. Therefore,

$\sum_{J}\Bigl{(}1+\sum_{b_{2}>C^{\prime}b_{1}}u_{N,\lambda,I,J}(b_{2})\Bigr{)}\ll\Bigl{(}\frac{\frac{hN^{2}}{MB^{3}}}{\frac{N}{M^{1/3}B^{1/3}}}+1\Bigr{)}\frac{N/L}{\frac{cN^{3}}{MB^{3}}}\ll\frac{hM^{1/3}B^{1/3}}{LN}+\frac{MB^{3}}{LN^{2}}$

and, summing over intervals $I$ in $(B,2B]$ and over dyadic intervals $(B,2B]$
up to $N\log^{1/2}x/h$ and combining with (45), we have

$\displaystyle\sum_{R<b_{1}\leq N\log^{1.2}x/h}$ $\displaystyle
t_{N,\lambda}(b_{1})=\sum_{B}\sum_{b_{1}\in(B,2B]}t_{N,\lambda}(b_{1})$
$\displaystyle\ll\sum_{B}\Bigl{(}\frac{hM^{1/3}B^{1/3}}{LN}+\frac{MB^{3}}{LN^{2}}\Bigr{)}\frac{B}{\frac{0.0001N}{M}}+\frac{x^{1/4}}{LN^{1/4}}$
$\displaystyle\ll\frac{M^{4/3}\log^{1.6}x}{LN^{2/3}h^{1/3}}+\frac{M^{2}N\log^{4.8}x}{Lh^{4}}+\frac{x^{1/4}}{LN^{1/4}}=\frac{x^{4/3}\log^{1.6}x}{LN^{10/3}h^{1/3}}+\frac{x^{2}\log^{4.8}x}{LN^{3}h^{4}}+\frac{x^{1/4}}{LN^{1/4}}\leq\frac{h}{13}$

as $x^{1/6}<h<x^{1/5}$ and $N\geq x^{2/3}\log^{2}x/h^{5/3}$ by choosing $L$
sufficiently large. Together with (4), we have the proposition.

##  8 Proof of Proposition 4

From this point onward, we will consider quadruples $(m+a_{1})(n-b_{1})^{2}$,
$mn^{2}$, $(m-a_{2})(n+b_{2})^{2}$ and $(m-a_{2}-a_{3})(n+b_{2}+b_{1})^{2}$
with $n-b_{1},n+b_{2}\in T_{N,\lambda}(b_{1})$ for some integers $b_{1}>R$ and
$b_{2}>0$ (i.e. $b_{1}=b_{3}$). If $b_{2}\leq C^{\prime}b_{1}$ for some large
constant $C^{\prime}>0$, then (6) and (33) imply

$1\leq D^{\prime}\ll\frac{Mb_{1}^{4}}{N^{3}}\;\;\text{ or
}\;\;b_{1}\gg\frac{N^{3/4}}{M^{1/4}},$

and the contribution from such quadruples with $b_{2}\leq C^{\prime}b_{1}$ is

$\leq\sum_{b_{1}\gg
N^{3/4}/M^{1/4}}t_{N,\lambda}(b_{1})\ll\frac{N/L}{N^{3/4}/M^{1/4}}=\frac{x^{1/4}}{LN^{1/4}}.$
(46)

Now, we focus on those quadruples with $b_{2}>C^{\prime}b_{1}$. We break
$(N,\lambda N]$ into subintervals of length $\frac{cN^{3}}{Mb_{1}^{3}}$ for
some sufficiently small $c>0$ so that $D_{4}=0$ by (6). Then (6) further
implies

$\frac{Mb_{1}^{3}b_{2}}{N^{3}}\ll\frac{h}{N}\;\;\text{ or
}\;\;b_{2}\ll\frac{hN^{2}}{Mb_{1}^{3}}$

which gives an upper bound of the distance between any two members in
$T_{N,\lambda}(b_{1})$ when $n$ is in such a subinterval. Meanwhile, by (33),
$D^{\prime}\geq 1$ which implies

$\frac{Mb_{1}b_{2}^{3}}{N^{3}}\gg 1\;\;\text{ or
}\;\;b_{2}\gg\frac{N}{M^{1/3}b_{1}^{1/3}},$

a lower bound for such distance. If $b_{1}\leq N^{1/5-\delta}$, then the
contribution

$\sum_{b_{1}\leq N^{1/5-\delta}}t_{N,\lambda}(b_{1})\ll\sum_{b_{1}\leq
N^{1/5-\delta}}\Bigl{(}\frac{\frac{hN^{2}}{Mb_{1}^{3}}}{\frac{N}{M^{1/3}b_{1}^{1/3}}}+1\Bigr{)}\frac{N/L}{\frac{cN^{3}}{Mb_{1}^{3}}}\ll\frac{hM^{1/3}}{N^{11/15+4\delta/3}}+\frac{M}{N^{6/5+4\delta}}\ll\frac{hx^{1/3}}{N^{7/5+4\delta/3}}+\frac{x}{N^{16/5+4\delta}}.$
(47)

So, from now on, we only consider $b_{1}>Q=:\max(N^{1/5-\delta},R)$. Suppose
$N\leq C_{\alpha}x^{\alpha}$ for some $\theta<\alpha<2/5$ with some constant
$C_{\alpha}>0$. For a fixed $b_{1}>Q$, define

$U_{N,\lambda,b_{1}}(b_{2}):=\\{n:n\text{ and }n+b_{2}\text{ are consecutive
elements in }T_{N,\lambda}(b_{1})\\},\text{ and
}u_{N,\lambda,b_{1}}(b_{2}):=|U_{N,\lambda,b_{1}}(b_{2})|.$

Clearly, one has

$\displaystyle\sum_{Q<b_{1}\leq
N\log^{1.2}x/h}t_{N,\lambda}(b_{1})\leq\sum_{Q<b_{1}\leq
N\log^{1/2}x/h}\Bigl{(}1+\sum_{b_{2}}u_{N,\lambda,b_{1}}(b_{2})\Bigr{)}$
$\displaystyle\ll$
$\displaystyle\frac{N}{h}+\frac{x^{1/4}}{LN^{1/4}}+\frac{hx^{1/3}}{N^{7/5+4\delta/3}}+\frac{x}{N^{16/5+4\delta}}+\sum_{Q<b_{1}\leq
N\log^{1.2}x/h}\sum_{b_{2}>C^{\prime}b_{1}}u_{N,\lambda,b_{1}}(b_{2})$ (48)

by (46) and (47). We break the double sum in (8) to the following two pieces:

$\sum_{Q<b_{1}\leq
N\log^{1.2}x/h}\mathop{\sum_{b_{2}>C^{\prime}b_{1}}}_{b_{2}>N^{\beta}/b_{1}}u_{N,\lambda,b_{1}}(b_{2})+\sum_{Q<b_{1}\leq
N\log^{1.2}x/h}\mathop{\sum_{b_{2}>C^{\prime}b_{1}}}_{b_{2}\leq
N^{\beta}/b_{1}}u_{N,\lambda,b_{1}}(b_{2}):=S_{1}+S_{2}.$ (49)

In $S_{1}$, $b_{2}>C^{\prime}b_{1}$ and $b_{2}>N^{\beta}/b_{1}$. Again, we
break $(N,\lambda N]$ into subintervals of length $\frac{cN^{3}}{Mb_{1}^{3}}$
for some sufficiently small $c>0$ so that $D_{4}=0$ by (6). Then (6) further
implies

$\frac{Mb_{1}^{2}}{N^{3-\beta}}\ll\frac{Mb_{1}^{3}b_{2}}{N^{3}}\ll\frac{h}{N}$
(50)

which gives

$b_{1}\leq\frac{C_{1}h^{1/2}N^{1-\beta/2}}{M^{1/2}}$ (51)

for some $C_{1}>0$. As $b_{2}>C^{\prime}b_{1}$, (6) and (21) give

$\frac{Mb_{1}b_{2}^{3}}{N^{3}}\gg D^{\prime}\geq 1\;\;\text{ and
}\;\;b_{2}\gg\frac{N}{M^{1/3}b_{1}^{1/3}}$ (52)

for sufficiently large $C^{\prime}>0$. Moreover, (50) implies

$b_{2}\ll\frac{hN^{2}}{Mb_{1}^{3}}.$ (53)

Therefore, combining (8), (51), (50), (52) and (53), we obtain

$\displaystyle S_{1}\ll$ $\displaystyle\sum_{Q<b_{1}\leq
C_{1}h^{1/2}N^{1-\beta/2}/M^{1/2}}\Bigl{(}\frac{\frac{hN^{2}}{Mb_{1}^{3}}}{\frac{N}{M^{1/3}b_{1}^{1/3}}}+1\Bigr{)}\frac{N/L}{\frac{cN^{3}}{Mb_{1}^{3}}}\ll\frac{h^{5/3}}{LM^{1/3}N^{(2\beta-1)/3}}+\frac{h^{2}N^{2-2\beta}}{LM}$
$\displaystyle\ll$
$\displaystyle\frac{hN^{1-2\beta/3}}{Lx^{1/3-2\theta/3}}+\frac{hN^{4-2\beta}}{Lx^{1-\theta}}\ll\frac{hN^{4-2\beta}}{Lx^{1-\theta}}$
(54)

as one can check that
$\frac{N^{1-2\beta/3}}{x^{1/3-2\theta/3}}\leq\frac{N^{4-2\beta}}{x^{1-\theta}}$
when $N\leq C_{\alpha}x^{\alpha}$ and $\alpha<2/5$. By (6) with $b_{1}=b_{3}$,
we have $1\leq a_{1}-a_{3}\ll\frac{mb_{1}b_{2}}{n^{2}}$. If $b_{1}b_{2}\leq
N^{\beta}$, this implies

$\frac{N^{4}}{x}\ll\frac{N^{2}}{M}\ll b_{1}b_{2}\leq N^{\beta}\;\;\text{ or
}\;\;N\ll x^{\frac{1}{4-\beta}}$

In view of the first lower bound for $N$ in (7), this means that the sum
$S_{2}$ is empty if $C_{\alpha}^{\prime}$ is sufficiently large, and we have
the proposition by (8), (49) and (8).

One can modify the above argument to account for the other bounds in (7). By
(6) and (6) with $b_{1}=b_{3}$, we have

$\frac{N}{b_{2}^{2}}\leq\frac{D^{\prime}N}{b_{2}^{2}}\ll
a_{1}-a_{3}\ll\frac{Mb_{1}b_{2}}{N^{2}}\leq\frac{M}{N^{2-\beta}}$

when $b_{1}b_{2}\leq N^{\beta}$. This implies

$\frac{N^{(3-\beta)/2}}{M^{1/2}}\ll b_{2}\leq\frac{N^{\beta}}{b_{1}}\;\;\text{
or }\;\;b_{1}\ll\frac{M^{1/2}}{N^{3/2-3\beta/2}}.$

However, we have $b_{1}>Q\geq R$. This gives

$\frac{N^{2/3}}{M^{1/3}}\ll\frac{N^{1/2}}{M^{3/2-3\beta/2}}\;\;\text{ or
}\;\;N\ll x^{\frac{5}{23-9\beta}}.$ (55)

In view of the second lower bound for $N$ in (7), $S_{2}$ is empty if
$C_{\alpha}^{\prime}$ is sufficiently large which gives the proposition.

Similar to (55), one may argue with $b_{1}>Q\geq N^{1/5-\delta}$. Then

$N^{1/5-\delta}\ll\frac{M^{1/2}}{N^{3/2-3\beta/2}}\;\;\text{ or }\;\;N\ll
x^{\frac{1}{27/5-3\beta-2\delta}}.$

In view of the third lower bound for $N$ in (7), $S_{2}$ is empty if
$C_{\alpha}^{\prime}$ is sufficiently large which gives the proposition.

##  9 Proof of Proposition 2

Recall $\delta=0.0751$. Due to Propositions 1 and 3, we may consider
$x^{3/13}\asymp ch^{5/2}/x^{1/4}\leq N\leq x^{2/3}\log^{2}x/h^{2/3}\asymp
x^{9/26}\log^{2}x$ as $h=Cx^{5/26}$. Note that the first error term in (6) is
$\ll x^{0.16}$ as $N\ll x^{9/26}\log^{2}x$, the third error term in (6) is
$\ll h/x^{0.01}$ as $x^{3/13}\ll N$, and the fourth error term in (6) is $\ll
h/x^{0.00009}$ as $\delta=0.0751$ and $x^{3/13}\ll N$. We apply (4) (similarly
for below) and Proposition 4 (using first lower bound in (7)) with
$\alpha=0.34616$ and $\beta=0.8334$, and obtain

$S(C_{1}x^{0.3158},x^{2/3}/h^{5/3})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.3158})^{1/4}}+\frac{h}{Lx^{0.00004}}\ll\frac{hL}{\log^{0.2}x}$
(56)

for some large $C_{1}>0$ after adding contributions from $\lambda$-adic
intervals $(N,\lambda N]$. Next, we apply Proposition 4 (using second lower
bound in (7)) with $\alpha=0.3158$ and $\beta=0.7212$, and obtain

$S(C_{2}x^{0.3029},C_{1}x^{0.3158})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.3029})^{1/4}}+\frac{h}{Lx^{0.000002}}\ll\frac{hL}{\log^{0.2}x}$
(57)

for some large $C_{2}>0$. Next, we apply Proposition 4 (using second lower
bound in (7)) with $\alpha=0.3029$ and $\beta=0.6668$, and obtain

$S(C_{3}x^{0.2942},C_{2}x^{0.3029})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2942})^{1/4}}+\frac{h}{Lx^{0.00003}}\ll\frac{hL}{\log^{0.2}x}$
(58)

for some large $C_{3}>0$. Next, we apply Proposition 4 (using second lower
bound in (7)) with $\alpha=0.2942$ and $\beta=0.6274$, and obtain

$S(C_{4}x^{0.2882},C_{3}x^{0.2942})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2882})^{1/4}}+\frac{h}{Lx^{0.00005}}\ll\frac{hL}{\log^{0.2}x}$
(59)

for some large $C_{4}>0$. Next, we apply Proposition 4 (using second lower
bound in (7)) with $\alpha=0.2882$ and $\beta=0.5988$, and obtain

$S(C_{5}x^{0.2840},C_{4}x^{0.2882})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2840})^{1/4}}+\frac{h}{Lx^{0.00004}}\ll\frac{hL}{\log^{0.2}x}$
(60)

for some large $C_{5}>0$. Next, we apply Proposition 4 (using second lower
bound in (7)) with $\alpha=0.2840$ and $\beta=0.5781$, and obtain

$S(C_{6}x^{0.2810},C_{5}x^{0.2840})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2810})^{1/4}}+\frac{h}{Lx^{0.00005}}\ll\frac{hL}{\log^{0.2}x}$
(61)

for some large $C_{6}>0$. Next, we apply Proposition 4 (using second lower
bound in (7)) with $\alpha=0.2810$ and $\beta=0.5629$, and obtain

$S(C_{7}x^{0.2789},C_{6}x^{0.2810})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2789})^{1/4}}+\frac{h}{Lx^{0.00004}}\ll\frac{hL}{\log^{0.2}x}$
(62)

for some large $C_{7}>0$. Next, we apply Proposition 4 (using second lower
bound in (7)) with $\alpha=0.2789$ and $\beta=0.5521$, and obtain

$S(C_{8}x^{0.2773},C_{7}x^{0.2789})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2773})^{1/4}}+\frac{h}{Lx^{0.00005}}\ll\frac{hL}{\log^{0.2}x}$
(63)

for some large $C_{8}>0$. Next, we apply Proposition 4 (using second lower
bound in (7)) with $\alpha=0.2773$ and $\beta=0.5437$, and obtain

$S(C_{9}x^{0.2762},C_{8}x^{0.2773})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2762})^{1/4}}+\frac{h}{Lx^{0.00002}}\ll\frac{hL}{\log^{0.2}x}$
(64)

for some large $C_{9}>0$. Next, we apply Proposition 4 (using third lower
bound in (7)) with $\alpha=0.2762$ and $\beta=0.5379$, and obtain

$S(C_{10}x^{0.2751},C_{9}x^{0.2762})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2751})^{1/4}}+\frac{h}{Lx^{0.00002}}\ll\frac{hL}{\log^{0.2}x}$
(65)

for some large $C_{10}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2751$ and $\beta=0.5321$, and obtain

$S(C_{11}x^{0.2738},C_{10}x^{0.2751})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2738})^{1/4}}+\frac{h}{Lx^{0.00005}}\ll\frac{hL}{\log^{0.2}x}$
(66)

for some large $C_{11}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2738$ and $\beta=0.5251$, and obtain

$S(C_{12}x^{0.2722},C_{11}x^{0.2738})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2722})^{1/4}}+\frac{h}{Lx^{0.00003}}\ll\frac{hL}{\log^{0.2}x}$
(67)

for some large $C_{12}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2722$ and $\beta=0.5164$, and obtain

$S(C_{13}x^{0.2703},C_{12}x^{0.2722})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2703})^{1/4}}+\frac{h}{Lx^{0.00002}}\ll\frac{hL}{\log^{0.2}x}$
(68)

for some large $C_{13}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2703$ and $\beta=0.5060$, and obtain

$S(C_{14}x^{0.2680},C_{13}x^{0.2703})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2680})^{1/4}}+\frac{h}{Lx^{0.00003}}\ll\frac{hL}{\log^{0.2}x}$
(69)

for some large $C_{14}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2680$ and $\beta=0.4932$, and obtain

$S(C_{15}x^{0.2653},C_{14}x^{0.2680})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2653})^{1/4}}+\frac{h}{Lx^{0.00004}}\ll\frac{hL}{\log^{0.2}x}$
(70)

for some large $C_{15}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2653$ and $\beta=0.4778$, and obtain

$S(C_{16}x^{0.2621},C_{15}x^{0.2653})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2621})^{1/4}}+\frac{h}{Lx^{0.00001}}\ll\frac{hL}{\log^{0.2}x}$
(71)

for some large $C_{16}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2621$ and $\beta=0.4592$, and obtain

$S(C_{17}x^{0.2583},C_{16}x^{0.2621})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2583})^{1/4}}+\frac{h}{Lx^{0.000004}}\ll\frac{hL}{\log^{0.2}x}$
(72)

for some large $C_{17}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2583$ and $\beta=0.4366$, and obtain

$S(C_{18}x^{0.2539},C_{17}x^{0.2583})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2539})^{1/4}}+\frac{h}{Lx^{0.00003}}\ll\frac{hL}{\log^{0.2}x}+$
(73)

for some large $C_{18}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2539$ and $\beta=0.4095$, and obtain

$S(C_{19}x^{0.2487},C_{18}x^{0.2539})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2487})^{1/4}}+\frac{h}{Lx^{0.00003}}\ll\frac{hL}{\log^{0.2}x}$
(74)

for some large $C_{19}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2487$ and $\beta=0.3762$, and obtain

$S(C_{20}x^{0.2427},C_{19}x^{0.2487})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2427})^{1/4}}+\frac{h}{Lx^{0.00001}}\ll\frac{hL}{\log^{0.2}x}$
(75)

for some large $C_{20}>0$. Next, we apply Proposition 4 (using the third lower
bound in (7)) with $\alpha=0.2427$ and $\beta=0.3361$, and obtain

$S(C_{21}x^{0.2358},C_{20}x^{0.2427})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(x^{0.2358})^{1/4}}+\frac{h}{Lx^{0.00003}}\ll\frac{hL}{\log^{0.2}x}$
(76)

for some large $C_{21}>0$. Finally, we apply Proposition 4 (using the third
lower bound in (7)) with $\alpha=0.2358$ and $\beta=0.2874$, and obtain

$S(ch^{5/2}/x^{1/4},C_{21}x^{0.2358})\ll\frac{hL}{\log^{0.2}x}+\frac{h}{x^{0.00009}}+\frac{x^{1/4}}{L(ch^{5/2}/x^{1/4})^{1/4}}+\frac{h}{Lx^{0.00003}}\ll\frac{h}{Lc^{1/4}}.$
(77)

Combining (56) - (77), we have the proposition by picking $L$ sufficiently
large.

## References

  * [1] M. Filaseta, An elementary approach to short interval results for $k$-free numbers, J. Number Theory 30 (1988), 208–225.
  * [2] M. Filaseta, Short interval results for squarefree numbers, J. Number Theory 35 (1990), 128–149.
  * [3] M. Filaseta and O. Trifonov, On gaps between squarefree numbers, Analytic Number Theory, Proceedings of a Conference in Honor of Paul T. Bateman (Progress in Mathematics Series, Vol. 85), edited by Berndt, Diamond, Halberstam, and Hildebrand, Birkhäuser, Boston (1990), 235–253.
  * [4] M. Filaseta and O. Trifonov, On gaps between squarefree numbers II, J. London Math. Soc. (2) 45 (1992), 323–333.
  * [5] M. Filaseta and O. Trifonov, The distribution of fractional parts with applications to gap results in number theory, Proc. London Math. Soc. (2) 3 (1996), 241–278.
  * [6] E. Fogels, On the average values of arithmetic functions, Proc. Cambridge Philos. Soc. 37 (1941), 358–372.
  * [7] S.W. Graham and G. Kolesnik, On the difference between consecutive squarefree integers, Acta Arith. 49 (1987), 435–447.
  * [8] R.A. Rankin, Van der Corput’s method and the theory of exponent pairs, Quart. J. Math. Oxford 6 (1955), 147–153.
  * [9] H.E. Richert, On the difference between consecutive square-free numbers, J. London Math. Soc. 29 (1954), 16–20.
  * [10] K.F. Roth, On the gaps between squarefree numbers, J. London Math. Soc. 26 (1951), 263–268.
  * [11] P.G. Schmidt, Abschätzungen bei unsymmetrischen Gitterpunktproblemen, Dissertation zur Erlangung des Doktorgrades der Mathematisch-Naturwissenschaftlichen Fakultät der Gerog-August-Universität zu Göttingen, 1964.
  * [12] O. Trifonov, On the squarefree problem, C. R. Acad. Bulgare Sci. 41 (1988), 37–40.
  * [13] O. Trifonov, On the squarefree problem II, Mathematica Balcanika 3 (1989), 284–295.

Mathematics Department
Kennesaw State University
Marietta, GA 30060


# Geodesic Regression Characterizes
3D Shape Changes in the Female Brain During Menstruation

Adele Myers
UC Santa Barbara
<EMAIL_ADDRESS>Caitlin Taylor
UC Santa Barbara
<EMAIL_ADDRESS>Emily Jacobs
UC Santa Barbara
<EMAIL_ADDRESS>Nina Miolane
UC Santa Barbara
<EMAIL_ADDRESS>
###### Abstract

Women are at higher risk of Alzheimer’s and other neurological diseases after
menopause, and yet research connecting female brain health to sex hormone
fluctuations is limited. We seek to investigate this connection by developing
tools that quantify 3D shape changes that occur in the brain during sex
hormone fluctuations. Geodesic regression on the space of 3D discrete surfaces
offers a principled way to characterize the evolution of a brain’s shape.
However, in its current form, this approach is too computationally expensive
for practical use. In this paper, we propose approximation schemes that
accelerate geodesic regression on shape spaces of 3D discrete surfaces. We
also provide rules of thumb for when each approximation can be used. We test
our approach on synthetic data to quantify the speed-accuracy trade-off of
these approximations and show that practitioners can expect very significant
speed-up while only sacrificing little accuracy. Finally, we apply the method
to real brain shape data and produce the first characterization of how the
female hippocampus changes shape during the menstrual cycle as a function of
progesterone: a characterization made (practically) possible by our
approximation schemes. Our work paves the way for comprehensive, practical
shape analyses in the fields of bio-medicine and computer vision. Our
implementation is publicly available on GitHub.

##  1 Introduction

Women are more likely to experience Alzheimer’s, cognitive decline, and
navigational issues after menopause [2, 20]. Yet, topics relevant to female
brain health such as menstruation, pregnancy, menopause, and their associated
female sex hormone fluctuations only account for 0.3% of the neuroimaging
literature between 1995 and 2022 [17]. The hippocampal formation (a brain
structure) is an excellent diagnostic tool for investigating the connection
between female brain health and sex hormones, as it is the first cortical
region to harbor neuropathology in the progression to Alzheimer’s (causing
characteristic shape changes visible on magnetic resonance images (MRI)) [5,
13], and it is also very sensitive to sex hormone fluctuations [7]. Sex
hormone fluctuations significantly influence brain anatomy and function in
healthy subjects [18, 16, 15]. Enhancing our understanding of how hormonal
fluctuations affect the healthy brain is crucial to explaining why females are
later more at risk for neurological conditions after menopause [4]. We seek to
close this knowledge gap by starting with one question: How does the
hippocampal formation respond to monthly fluctuations in ovarian hormones
during the menstrual cycle?

Figure 1: During the menstrual cycle, the ovaries cyclically release female
sex hormones such as progesterone into the blood. We propose practical shape
analysis tools that quantify 3D shape changes in the brain that occur during
progesterone fluctuations, with a focus on the hippocampus: a structure
involved in memory and navigation. Visualization created from data by [18,
15].

Recent research has shown that certain substructures of the hippocampal
formation change their volume over the course of the menstrual cycle in
response to progesterone, but no significant volumetric change was found on
the whole-formation level [18]. Our 3D visualizations (Fig. 1) show that the
hippocampal formation does change its shape on a whole-formation level, but no
team has quantified how this effect depends on progesterone levels. This is
not surprising because quantifying 3D shape changes is technically
challenging, and is in fact is an active research area in itself in
mathematics and computer vision. For example, quantifying surface shape
changes as a function of a continuous variable like progesterone can
theoretically be performed through geodesic regression on Riemannian manifolds
[6, 19]: an extension of linear regression dedicated to shape spaces. However,
in its current form, geodesic regression on surface spaces is too slow for
practical use. Here, we bridge the gap between computer vision and clinical
neuroimaging by presenting a new practical method, a hybrid between geodesic
and linear regression, that allows us to quantify how the shape of the
hippocampal formation changes in response to progesterone.

#### Contributions

We offer several contributions that span the fields of machine learning,
differential geometry and clinical neuroimaging. First, in machine learning,
we introduce our hybrid geodesic-linear regression method: a faster geodesic
regression model which uses linear residuals instead of geodesic residuals in
its loss function and also uses a linear regression result to initialize its
geodesic regression optimization. Then, we perform extensive synthetic
experiments to offer rules of thumb for deciding between linear, geodesic, and
geodesic-linear regression. In differential geometry, these results give novel
intuition about the curvature of the nonlinear data space of surface shapes.
In clinical neuroimaging, these rules of thumb provide practitioners with
guidelines to decide on the speed-accuracy trade-off between the regression
types, revealing whether a surface mesh sequence can be adequately
characterized by the considerably faster linear or geodesic-linear regressions
without sacrificing accuracy. Finally, we apply our paradigm to real brain
magnetic resonance images (MRIs) of a female brain through the menstrual
cycle. We characterize, for the first time, shape changes in the female
hippocampal formation as a function of progesterone.

##  2 Related Works

Consider a series of hippocampal shapes, with the surface of each shape
described as a mesh extracted via segmentation from a full-brain MRI (Fig. 1).
The shapes of these discrete surfaces (meshes) can be described either
extrinsically or intrinsically. The extrinsic approach with the Large
Deformation Diffeomorphic Metric Mapping framework [3] represents a surface
shape in 3D space by the amount of deformation that one needs to apply on the
whole 3D ambient grid containing the surface in order to deform a reference
shape into the surface shape of interest. By contrast, the intrinsic approach
only deforms the surface itself [1], and hence provides us with two
advantages: (i) intuitive deformations that can be discussed with
neuroscientists, (ii) higher computational efficiency with up to 10x
acceleration [1]. We focus here on the intrinsic approach.

#### Analysis of Parameterized Surfaces

In the intrinsic approach, surfaces can be either parameterized or
unparameterized. Each mesh in a dataset of parameterized surfaces is
constrained to have a consistent structure: the vertices of meshes in the
dataset have one-to-one correspondences. By contrast, datasets of
unparameterized surfaces relax this constraint such that statistical analyses
can be performed independently of the number and indexation of the mesh
vertices. While theoretically grounded, the computational complexity of this
approach makes it unpractical for geodesic regression. As an example,
statistical analyses within this framework are often limited to population
average computations and machine learning algorithms that rely only on
distances such as multidimensional scaling or k-means clustering [1, 9, 12].

Therefore, we consider the scenario where each mesh in the dataset is first
parameterized to match a reference parameterization, and then statistical
analysis such as regression is performed. In this scenario, the first natural
choice is to consider linear regression on the set of 3D coordinates of the
vertices. Linear regression has the advantage of being conceptually simple
while enjoying analytical solutions: it will be our first baseline. However,
it has the drawback that it does not enjoy parameterization-invariance. In
other words, the distance (or dissimilarity) between two surfaces may change
if we choose another reference parameterization, which may change the
regression results. In the context of clinical application where the ultimate
goal is human health, it is thought that we cannot afford such inconsistency.
The alternative is to consider parameterization-invariant distances between
parameterized surfaces. In differential geometry, this can be achieved by
equipping the space of parameterized surfaces with a Riemannian metric that is
invariant under reparameterizations [10, 9]. This process however, turns the
data space into a nonlinear manifold, where linear regression needs to be
generalized to geodesic regression [6, 19]. Geodesic regression does not enjoy
a closed-form estimator and is typically computationally expensive. It will be
our second baseline.

Figure 2: A surface is represented by a function $q:M\to\mathbb{R}^{3}$ that
maps parameters $(u,v)\in M$ to points in 3D space $q(u,v)\in\mathbb{R}^{3}$
(top row). Its parameterization can be changed by applying a diffeomorphism
$\phi$ to the domain $M$ before mapping to $\mathbb{R}^{3}$ (bottom row).

#### Computational Challenges of Geodesic Regression

Geodesic regression [6, 19] for parameterization-invariant Riemannian metrics
presents unique computational challenges. In the Riemannian framework,
calculating a single geodesic requires a computationally expensive approach.
Geodesic regression solves an optimization problem by minimizing a mean square
error (MSE) loss function. The MSE requires the computations of $n+1$
geodesics at each iteration: 1 geodesic representing the generative model, and
$n$ geodesics required to compute the residuals in the MSE, where $n$ is the
number of surfaces (meshes) in the dataset. We observe that works developing
the Riemannian framework limit the number of geodesic computations required
for their analysis, and do not perform any form of regression. Kilian et al.
[11] focus on geodesic interpolation or extrapolation of parameterized
surfaces and do not study regression. Kurtek et al. [12] limit their
experimental analysis to computing geodesics between pairs of unparameterized
surfaces and performing clustering. Jermyn et al. [10] compute geodesics
between pairs of parameterized surfaces and provide a classification
experiment. Hartman et al [9] estimate population averages and perform
dimension reduction with multidimensional scaling (MDS) and tangent PCA for
parameterized and unparameterized surfaces. Bauer et al [1] compute geodesics
and the population averages of unparameterized surfaces, together with multi-
dimensional scaling and k-means clustering. We suspect that the authors did
not perform geodesic regression in these works because of their computational
costs, which we investigate here.

##  3 Background

This section presents the mathematical background necessary to formulate our
approximation schemes for geodesic regression on the shape space of
(hippocampal) surfaces. We refer to [9, 8] for additional details.

#### A. Riemannian Metrics and Geodesics

We first introduce concepts in differential geometry necessary for geodesic
regression. A Riemannian metric on a smooth manifold $\mathcal{N}$ is a family
$\left(G_{p}\right)_{p\in\mathcal{N}}$ of inner products on each tangent space
$T_{p}\mathcal{N}$, such that $G_{p}$ depends smoothly on the point
$p\in\mathcal{N}$. Any Riemannian metric $G$ yields a notion of distance
between points $q_{0},q_{1}$ on $\mathcal{N}$. Specifically, if
$\gamma:[0,1]\rightarrow\mathcal{N}$ is a smooth trajectory on $\mathcal{N}$
with velocity vector at $t\in[0,1]$ denoted as $\dot{\gamma}_{t}\in
T_{\gamma(t)}\mathcal{N}$, its length is defined as
$L_{\gamma}=\int_{0}^{1}\sqrt{G(\dot{\gamma}_{t},\dot{\gamma}_{t})_{\gamma_{t}}}dt$
and the distance between any two points $q_{0},q_{1}\in\mathcal{N}$ is given
by $d(q_{0},q_{1})=\inf_{\gamma:\gamma(0)=q_{0},\gamma(1)=q_{1}}L_{\gamma}$.

A geodesic between two points $q_{0},q_{1}$ that are “close” in $\mathcal{N}$
is defined as a trajectory $\gamma$ that locally realizes the shortest
distance between $q_{0}$ and $q_{1}$. Intuitively, a geodesic is the
generalization to manifolds of the concept of a straight line in vector
spaces. While some manifolds enjoy analytical expression for their geodesics,
this is not case for the manifold of (hippocampal) surface shapes that we will
consider here. Thus, geodesics will need to be computed numerically.

To this aim, geodesics are expressed as the solutions of the geodesic
equation, which is an ordinary differential equation (ODE) which can be
written in local coordinates as:

$\ddot{\gamma}^{k}(t)+\Gamma_{ij}^{k}\dot{\gamma}^{i}(t)\dot{\gamma}^{j}(t)=0,$
(1)

for all times $t\in[0,1]$ where $\Gamma_{ij}^{k}$ are the Christoffel symbols
associated with the Riemannian metric. Solving this ODE provides numerical
solutions for geodesics.

To perform geodesic regression, we will also need two additional operations,
called Exp and Log, which we define here. The map
$(q,v)\mapsto\gamma_{q,v}(1)$ defined for $(q,v)\in\mathcal{I}\times
T_{q}\mathcal{I}$ is called the exponential map (Exp) and essentially computes
the point $\gamma_{q,v}(1)$ after following the geodesic of initial point
$q\in\mathcal{N}$ and initial velocity $v\in T_{q}\mathcal{N}$. The inverse of
the Exp map on its injectivity domain is called the logarithm map (Log).

#### B. Surfaces and Their Parameterizations

A continuous surface can be described by a function $q:M\to\mathbb{R}^{3}$,
where $M$ is a two-dimensional space of parameters $(u,v)\in M$ that
parameterize the 3D points $q(u,v)\in\mathbb{R}^{3}$ on the surface.
Intuitively, the function $q$ deforms the space of parameters $M$ to give the
surface its distinct shape, e.g., the ellipsoid shown in the top row of Fig.
2. Mathematically, $q$ is required to be an oriented smooth mapping in
$C^{\infty}(M,\mathbb{R}^{3})$ that is also regular in the sense that its
differential $dq$ is injective everywhere on $M$.

The parameterization of a surface refers to the placement of points on the
surface. If we define one surface as $q:M\to\mathbb{R}^{3}$, then we can
describe the same surface with a different parameterization by
$q\circ\phi:M\to\mathbb{R}^{3}$, where $\phi$ is an orientation-preserving
diffeomorphism of $M$. Intuitively, $\phi$ smoothly deforms the placement of
parameters on the domain $M$, which in turn smoothly changes the placement of
points in the co-domain, as shown with rainbow colors in the bottom row of
Fig. 2. The change of parameterization $\phi$ does not change the shape of the
surface, which is an ellipsoid in both rows of Fig. 2.

#### C. Space of Surfaces

The space of surfaces is denoted $\mathcal{I}\subset
C^{\infty}(M,\mathbb{R}^{3})$. The space $\mathcal{I}$ is an infinite
dimensional manifold immersed in the infinite dimensional vector space
$C^{\infty}(M,\mathbb{R}^{3})$. The Riemannian metric we choose to equip the
manifold $\mathcal{I}$ with defines the distance between its points
$q_{0},q_{1}\in\mathcal{I}$ and thus the notion of dissimilarity between the
two surfaces $q_{0},q_{1}$. We consider the second-order Sobolev metric [9]:

$\displaystyle G_{q}(h,k)$ $\displaystyle=\int_{M}\left(a_{0}\langle
h,k\rangle+a_{1}g_{q}^{-1}\left(dh_{m},dk_{m}\right)\right.$ (2)
$\displaystyle+b_{1}g_{q}^{-1}\left(dh_{+},dk_{+}\right)+c_{1}g_{q}^{-1}\left(dh_{\perp},dk_{\perp}\right)$
$\displaystyle\left.+d_{1}g_{q}^{-1}\left(dh_{0},dk_{0}\right)+a_{2}\left\langle\Delta_{q}h,\Delta_{q}k\right\rangle\right)\operatorname{vol}_{q},$

where $h,k$ are tangent vectors at point $q\in\mathcal{N}$; $g_{q}^{-1}$ is
the pullback metric from $\mathbb{R}^{3}$ that defines distances on the
surface $q$ itself; $\Delta_{q}$ is the Laplacian induced by $q$;
$dh_{m},dh_{+},dh_{\perp},dh_{0}$ are orthogonal vector-valued one-forms and
$\text{vol}_{q}$ is the surface area measure of $q$. The scalars
$a_{0},a_{1},a_{2},b_{1},c_{1},d_{1}$ are weighting parameters that define
distance between two surfaces based on how they are sheared, scaled, bent, or
parameterized with respect to each other.

The choice of second-order Sobolev metric is motivated by the following facts.
First, the zero-order and first-order Sobolev metrics yield less stable
results in geodesic interpolation between complex 3D shapes [9]. Second, the
weighting parameters $a_{0},a_{1},a_{2},b_{1},c_{1},d_{1}$ defining the
second-order Sobolev metric in Eq. (2) can be linked to observable physical
deformations (shearing, bending, etc) which helps with intuitively comparing
physical objects. Last, the metric in Eq. (2) yields a distance that is
rotation and reparameterization invariant [9]. In other words, if all the
surfaces in the dataset are rotated and reparameterized in the same way, i.e.,
using the same rotation matrix and reparameterization diffeomorphism $\phi$,
then their pairwise distances are unchanged. We note that this property is
practical only if we first assume that all the surfaces (are oriented and)
have valid point-to-point correspondences.

Figure 3: Distances in surface space vs shape space. $q_{1}$ and $q_{2}$ are
two surfaces with different parameterization and different shape. The distance
given by the second-order Sobolev metric 2 measures both the parameterization
and shape differences. The shape space distance 3 only measures difference in
shape.

#### D. Space of Surface Shapes

In the space of surfaces $\mathcal{I}$, if two surfaces have the same shape
but different orientations or parameterizations, they correspond to different
points. By contrast, we introduce the space of surface shapes where two
surfaces with the same shape correspond to the same point, regardless of
differences in their orientation or parameterization. Mathematically, the
space of surface shapes is defined as the quotient space:
$\mathcal{I}/(\text{Rot}(\mathbb{R}^{3})\times\text{Diff}(M))$ —see [9] for
details. For simplicity, we consider the case of parameterizations with the
shape space $\mathcal{S}=\mathcal{I}/\text{Diff}(M)$ while the case of
orientations can be treated similarly.

In the shape space $\mathcal{S}$, the distance between two surface shapes
$q_{1}$ and $q_{2}$ is given by:

$d^{\mathcal{S}}(q_{1},q_{2})=\text{inf}_{\phi}d(q_{1},q_{2}\circ\phi)\\\
=d(q_{1},q_{2}^{\prime}),$ (3)

where $\phi$ represents a choice in parameterization. In Eq. (3), the
parameterization of $q_{2}$ is varied until the second-order Sobolev distance
$d$ in Eq. (2) between $q_{1}$ and $q_{2}$ reaches an infimum as shown in Fig.
3. This operation matches the parameterization of $q_{2}$ to the
parameterization of $q_{1}$ so that any remaining discrepancy between them is
due to difference in shape, rather than difference in parameterization.
Ideally, we would perform our geodesic regression methods directly in the
shape space $\mathcal{S}$. However, the high computational cost of this
approach leads us to instead compute in the surface space $\mathcal{I}$ after
choosing a reference parameterization that corresponds to the first
hippocampal surface of our dataset.

##  4 Methods

Figure 4: Overview: Approximation schemes for geodesic regression
approximation. $\delta$-test: if the residual magnitudes are small compared to
the curvature of the manifold, we can use geodesic regression with linear
residuals (GRLR). $\Delta$-test: if the distance covered by the data set is
small compared to the curvature, we use linear regression (LR). Reducing
geodesic regression (GR) to either GRLR or LR provides up to four orders of
magnitude speed-up, while sacrificing little accuracy.

We seek to quantify the anatomical changes in the hippocampal formation that
emerge from progesterone variations during the menstrual cycle. To achieve
this, we propose approximations to geodesic regression on the space of 3D
brain shapes that make it computationally fast enough for practical use. We
further propose rules of thumb for determining when each approximation can be
used, as summarized in Fig. 4.

###  4.1 Linear Regression

#### Model

Linear regression (LR) models the relationship between an independent variable
$X\in\mathbb{R}$ and the dependent variable $Y$ taking values in
$\mathbb{R}^{D}$ as:

$Y=\alpha+X\beta+\epsilon,$ (4)

where $\alpha\in\mathbb{R}^{D}$ is the intercept, $\beta\in\mathbb{R}^{D}$ is
the slope, and $\epsilon$ represents the noise.

#### Loss

Given data $\left(x_{i},y_{i}\right)\in$ $\mathbb{R}\times\mathbb{R}^{D}$, for
$i=1,\ldots,n$, we fit the linear regression model through least squares,
i.e., we compute the estimates for the intercept and slope
$\hat{\alpha},\hat{\beta}$ as:

$(\hat{\alpha},\hat{\beta})=\arg\min_{(\alpha,\beta)}\frac{1}{2}\sum_{i=1}^{n}\left\|y_{i}-\hat{y}_{i}\right\|^{2}\text{
for }\hat{y}_{i}=\alpha+x_{i}\beta,$ (5)

which minimizes the summed squared magnitude of the (linear) residuals:
$y_{i}-\alpha-x_{i}\beta$, for $i=1,...,n$.

#### Learning

Importantly for computational purposes, this minimization has an analytical
solution given by the normal equations $\hat{\beta}=\frac{\frac{1}{n}\sum
x_{i}y_{i}-\bar{x}\bar{y}}{\sum x_{i}^{2}-\bar{x}^{2}}$ and
$\hat{\alpha}=\bar{y}-\bar{x}\hat{\beta}$, where $\bar{x}$ and $\bar{y}$ are
the sample means of the $x_{i}$ and $y_{i}$, respectively. We will use linear
regression as our first baseline, where $X$ is the level of progesterone, and
$Y$ is the hippocampal surface discretized as a mesh, which takes values in
$\mathbb{R}^{N\times 3}$ where $N$ is the number of mesh vertices.

###  4.2 Geodesic Regression

#### Model

Geodesic regression (GR) [6, 19] models the relationship between an
independent variable $X\in\mathbb{R}$ and the dependent variable $Y$, whose
values lie on a manifold $\mathcal{N}$, as:

$Y=\operatorname{Exp}(\operatorname{Exp}(p,Xv),\epsilon),$ (6)

where $\epsilon$ is noise in the tangent space at $\operatorname{Exp}(p,Xv)$,
and Exp is the operation defined in the previous section. Note that when the
manifold of interest is $\mathcal{N}=\mathbb{R}^{D}$, the exponential operator
simplifies to addition: $\operatorname{Exp}(p,v)=$ $p+v$. Consequently, the
geodesic regression generative model simplifies to the linear regression
generative model of Eq. (4) with $p=\alpha$ and $v=\beta$. We also note that
the exponential operation appears twice: to model the geodesic itself, and to
model the noise $\epsilon$. In what follows, we consider geodesic regression
on the manifold $\mathcal{N}=\mathcal{I}$ equipped with a second-order Sobolev
metric from Eq. (2).

#### Loss

Given data $\left(x_{i},y_{i}\right)\in\mathbb{R}\times\mathcal{I}$, for
$i=1,\ldots,n$, we seek to learn estimates of the intercept and slope
$(p,v)\in\mathcal{I}\times T_{p}\mathcal{I}$. In the manifold setting, the
loss function associated with the geodesic given by $(p,v)$ is minimized as:

$\displaystyle(\hat{p},\hat{v})$
$\displaystyle=\arg\min_{(p,v)}\frac{1}{2}\sum_{i=1}^{n}d\left(y_{i},\hat{y}_{i}\right)^{2},$
(7)
$\displaystyle=\arg\min_{(p,v)}\frac{1}{2}\sum_{i=1}^{n}\|\text{Log}(\hat{y}_{i},y_{i})\|_{\hat{y}_{i}}^{2},$
(8) $\displaystyle\text{for
}\hat{y}_{i}=\operatorname{Exp}\left(p,x_{i}v\right).$ (9)

We compute estimates for intercept and slope $(\hat{p},\hat{v})$ which
minimize the summed squared magnitude of the (geodesic) residuals
$\text{Log}(\operatorname{Exp}\left(p,x_{i}v\right),y_{i})$ for $i=1,...,n$.
The geodesic residuals differ from the linear residuals as they are calculated
with exponentials and logarithms instead of additions and subtractions.

#### Learning

In contrast to linear regression, the least squares problem of Eq. (7) above
does not have an analytical solution for general manifolds $\mathcal{I}$.
Instead, we need to compute the estimates of the intercept and slope with
gradient descent, which is typically computationally expensive. Gradient
descent comes in two flavors depending on the strategy used to compute the
gradient, which can be either a Riemannian gradient as originally proposed in
[6] or an extrinsic gradient. The Riemannian gradient writes [6]:

$\displaystyle\nabla_{p}l=-\sum_{i=1}^{N}d_{p}\operatorname{Exp}\left(p,x_{i}v\right)^{\dagger}\epsilon_{i},$
$\displaystyle\nabla_{v}l=-\sum_{i=1}^{N}x_{i}d_{v}\operatorname{Exp}\left(p,x_{i}v\right)^{\dagger}\epsilon_{i},$

where $l$ is the loss function,
$\epsilon_{i}=\text{Log}\left(\operatorname{Exp}\left(p,x_{i}v\right),y_{i}\right)$
are the residuals, $d_{v}$ and $d_{p}$ are derivatives, and $\dagger$ denotes
the adjoint. In the general case, the expression of these derivatives and
their respective adjoint operators are not known although they can be derived
analytically for some manifolds as in [19]. However, to the best of our
knowledge, no such formula exists for shape spaces of parameterized surfaces,
so we use a numerical approach.

###  4.3 Why is Geodesic Regression Slow?

The geodesic regression optimization is slow due to the Exp and Log maps in
Eq. (7). Computation of the exponential and logarithm maps do not enjoy an
analytical expression for the manifold that we are interested in, and neither
do their differentials. Consequently, we compute them only numerically, as
implemented in Geomstats [14] as follows.

For the exponential map, we consider the geodesic equation as a coupled system
of first-order ODEs:

$\left\\{\begin{array}[]{l}v(t)=\dot{\gamma}(t)\\\
\dot{v}(t)=f(v(t),\gamma(t),t)\end{array}\right.$

where $f$ is a smooth function given by Eq. (1) and the state variable is
$(\gamma(t),\dot{\gamma}(t))$. Given initial conditions, we use a first-order
forward Euler scheme to integrate this system. For a given step $dt$ we
compute:

$v(t+dt)=v(t)+\dot{v}(t)dt=v(t)+f(v(t),\gamma(t),t)dt.$ (10)

Introducing the parameter $n_{steps}$, if we integrate this geodesic equation
between $t=0$ and $t=1$ in $n_{\text{steps}}$ then we use
$dt=(n_{\text{steps}})^{-1}$. Consequently, the parameter $n_{\text{steps}}$
controls the numerical precision of the computation of the exponential map.
The computation of Exp is slow due to this numerical integration.

For the logarithm map, we solve the optimization problem in $v$:

$\min d^{2}\left(\operatorname{Exp}(p,v),q\right),$

that represents the fact that Log is the inverse map of Exp. This minimization
is solved by gradient descent (GD) until a convergence tolerance is reached.
It uses scipy for minimization method and computes the gradient of the
exponential map with automatic differentiation. The computation of Log is slow
due to this optimization process.

###  4.4 Approximations Schemes with Rules of Thumb

Curved spaces are locally linear. Thus, if a data set falls on a “small”
portion the shape space, addition and subtraction will offer excellent
approximations of exponentials and logarithms and avoid costly computations.
To speed up geodesic regression, we propose two approximation schemes shown in
Fig. 4: (i) linear regression, and (ii) geodesic regression with linear
residuals, where (ii) represents a novel approach for geometric machine
learning—which we describe in the next section.

We also propose rules of thumb to determine when each of these methods will
yield sufficiently accurate approximations of geodesic regression. For (i), we
propose the $\Delta$-Test, which explores when the magnitude of the geodesic
length of the data set ($\Delta$) is small at the scale of the curvature of
the manifold (Fig. 4 middle). For (ii), we propose the $\delta$-Test, which
explores when the magnitude of the noise ($\delta$) is small at the scale of
the curvature of the manifold (Fig. 4 right). In the experiments section, we
explore the curvature of the space of 3D discrete surfaces to give numerical
values to these guidelines.

###  4.5 Geodesic Regression with Linear Residuals

#### Model

We propose geodesic regression with linear residuals (GRLR) to model the
relationship between an independent variable $X\in\mathbb{R}$, the noise-free
dependent variable taking values in a manifold, and the (noisy) dependent
variable $Y$ taking values in $\mathbb{R}^{D}$. In other words, we propose the
following generative model:

$Y=\operatorname{Exp}(p,Xv)+\epsilon,$ (11)

where $\operatorname{Exp}(p,Xv)$ is the noise-free dependent variable and
$\epsilon$ is the noise. The noise-free dependent variable is constrained to
be a surface in $\mathcal{I}$, and thus $Y$’s dependency on $X$ is modelled
using the Exp operation. However, in practical applications, the data’s noise
may push the data off of $\mathcal{I}$. Thus, in addition to its computational
gain, this generative model acknowledges the fact that there is no reason for
the noise to be constrained on the manifold.

#### Loss

Given data $\left(x_{i},y_{i}\right)\in\mathbb{R}\times\mathbb{R}^{D}$, for
$i=1,\ldots,n$, we fit this regression model through least squares, i.e., we
compute the estimates for the intercept and slope as:

$(\hat{p},\hat{v})=\arg\min_{(p,v)}\frac{1}{2}\sum_{i=1}^{n}\|\hat{y}_{i}-y_{i}\|^{2}\text{
for }\hat{y}_{i}=\operatorname{Exp}\left(p,x_{i}v\right),$ (12)

where the squared geodesic distance of Eq. 7 has been replaced by the squared
Euclidean distance, but we keep the exponential map defining the geodesic.

#### Learning

Like geodesic regression, this least squares problem still requires gradient
descent. The gradient can be computed as a Riemannian gradient or as an
extrinsic gradient. The Riemannian gradient is given by:

$\displaystyle\nabla_{p}l=-\sum_{i=1}^{N}d_{p}\operatorname{Exp}\left(p,x_{i}v\right)^{\dagger}\epsilon_{i},$
$\displaystyle\nabla_{v}l=-\sum_{i=1}^{N}x_{i}d_{v}\operatorname{Exp}\left(p,x_{i}v\right)^{\dagger}\epsilon_{i},$

where $l$ is the loss function,
$\epsilon_{i}=y_{i}-\operatorname{Exp}\left(p,x_{i}v\right)$ are the
residuals, $d_{v}$ and $d_{p}$ are derivatives, and $\dagger$ is the adjoint.
By avoiding the computation of $n$ logarithms, the learning process enjoys a
significant speed-up. Even within the use of the extrinsic gradient, the loss
function avoids the computations of these logarithms and is thus accelerated.
We note that we still need to compute the Exp by numerical integrations and
their derivatives which we do by automatic differentiation. Our implementation
is publicly available on GitHub.

##  5 Experiments

We investigate the curvature of the space of surfaces to quantify which
approximation scheme should be used on which dataset. Guided by this analysis,
we approximate geodesic regression on 3D hippocampal surfaces, giving the
first characterization of hippocampal formation’s shape change as a function
of progesterone.

###  5.1 Curvature Estimation with $\delta-$ and $\Delta-$ Tests

#### Simulations

We perform experiments on synthetic meshes to provide rules-of-thumb that help
the practitioner decide when linear regression or geodesic regression with
linear residuals can be used with little loss in accuracy on the space of
discrete surfaces.

Specifically, our experiments explore cases when lines can be used to
approximate geodesics. First, we compute both a line and a geodesic between
two meshes $q^{\text{start}}$ and $q^{\text{end}}$. Then, we compare the
meshes along the line in $\mathbb{R}^{N\times 3}$ to the meshes along the
geodesic in $\mathcal{I}$, where $N$ is the number of vertices in the 3D
meshes. We have either $n=5$ or $n=10$ meshes along each sequence. The start
mesh $q^{\text{start}}$ is an ellipsoid whose principal axes have length 2, 2,
and 3. The end mesh is a deformed version $q^{\text{end}}$ of the reference
$q^{\text{start}}$, where the amount of deformation is controlled by a factor
that we vary in $\\{1\%,10\%,50\%,100\%\\}$. This deformation factor indicates
by how much each vertex in $q^{\text{end}}$ has been moved compared to the
vertex’s position in $q^{\text{start}}$ and is given as a percentage of the
diameter of $q^{\text{start}}$. $q^{end}$ is generated by adding isotropic
Gaussian noise to each vertex in $q^{start}$. In other words, the standard
deviation of the Gaussian noise is: $\sigma=\text{deformation}\times D$ where
$D$ is the diameter of the mesh.

We determine how much the geodesic and the line between $q^{\text{start}}$ and
$q^{\text{end}}$ differ by computing the root mean square deviation (RMSD)
between the geodesic mesh sequence and the line mesh sequence, which we then
normalize by the diameter $D$ of the mesh:

$\text{RMSD}=\frac{1}{D}\sqrt{\frac{1}{TN}\sum_{t=1}^{T}\sum_{j=1}^{N}\|v_{tj}^{\text{line}}-v_{tj}^{\text{geodesic}}\|^{2}},$
(13)

where $N$ is the number of vertices, $T$ the number of meshes in the sequence
(5 or 10) and $D$ the diameter. We also time the computation of the geodesic
and the line and report their ratio.

Figure 5: Accuracy-speed trade-off between a geodesic and its linear
approximation in joining two meshes $q^{\text{start}}$ and $q^{\text{end}}$.
The x-axis (error) shows how their meshes differ by computing distances
between their vertices. The y-axis (speed) shows the ratio of their
computational times. The color represents the deformation factor, i.e. how
deformed the mesh $q^{\text{end}}$ is from the reference mesh
$q^{\text{start}}$. The symbols represent the number of vertices in the
meshes. The number of steps is a parameter controlling the numerical
integration computing the geodesic. The results are similar when $n_{steps}$ =
5. Figure 6: Linear regression reveals hippocampal deformation associated with
an increase in progesterone during the menstrual cycle. The two rows show two
views of each 3D mesh along the sequence. For visualization purposes, a
coloring of the mesh is used to represent depth. While volumetric analyses did
not capture a volumetric change, our shape analysis reveals that an increase
in progesterone corresponds to a shear deformation of the hippocampal
formation.

####  $\delta$-Test

Fig. 5 shows that deformation factors of 1% yield errors below 0.05% of the
diameter of the shape (below 0.0005 on the figure). This is true across two
values of the number of steps used for the numerical computation of the
exponential map: $n_{steps}$ = 20, and $n_{steps}$ = 5 (see Section ).
Consequently, for our $\delta$-Test: when the measurement noise on the
vertices of the meshed shapes is expected to be less than 1% of the total mesh
diameter, we recommend using linear residuals instead of geodesic residuals.
In this case, we assess that the shape manifold can be approximated as linear
at the scale of residual length: the curvature is low compared to the
magnitude of the noise. By using linear residuals, we enjoy a considerable
speed-up, up to 14M $\times$ (for a number of steps of 20), as shown in Fig.
5, and up to 1.5M $\times$ (for a number of steps of 5) .

As a real-world example, consider $q^{\text{start}}$ as the mesh corresponding
to a true hippocampal shape at a given level of progesterone, and
$q^{\text{end}}$ as the mesh that we observe in practice after segmenting and
extracting the mesh from the MRI data. In this case, the “measurement noise”
is the MRI noise and the segmentation error. We expect measurement noise to
displace each vertex by around 1% of the total mesh diameter, since MRI images
have a good resolution and both segmentation and meshing algorithms are
reasonably accurate. Thus, brain MRIs are in the regime where geodesic
regression can enjoy considerable speed-ups by utilizing linear residuals.

####  $\Delta$-Test

Additionally, Fig. 5 shows that deformation factors of 10-50% yield an error
that is less than 10% of the diameter (below 0.1 on the figure). Consequently,
for our $\Delta$-Test: if the data set’s largest deformation between two
meshes is 10-50% of the diameter of the mesh, and if practitioners can
tolerate a maximum loss of accuracy of 10% in their results, then using linear
regression instead of geodesic regression allows them to significantly speed
up their pipeline. In this case, we assess that the shape manifold can be
approximated as linear on the scale of the data set: the curvature is low
compared to the spread of the data. Even more strikingly, the error decreases
from 10% to only 3% if practitioners consider meshes with hundreds of
vertices, as shown in Fig. 5.

###  5.2 Hippocampal Shape Change Characterization

#### Data set

We use a time-series of 3D brain images recorded from magnetic resonance
imaging (MRI): 11 images from 11 consecutive days, capturing the progesterone
peak of a single female subject’s natural menstrual cycle [15], as analyzed by
volumetric analyses in [18]. We choose to focus on the progesterone peak (11
days), as opposed to the full menstrual cycle (30 days) for simplicity. The
female also measured hormone levels in her blood in conjunction with each MRI
session.

#### Pre-processing

We align the 3D images to correct for the position and orientation of the
subject’s head in the MRI scanner, and to extract the surface of each
substructure of the hippocampal formation. The results of this pre-processing
are shown in Figure 1 where each sub-structure surface is color-coded and
shown at two different levels of progesterone (low and high). Here, we can
visually observe the hippocampal formation’s shape evolution, which our
analyses will seek to characterize with the slope and intercept learned from
the proposed approximation of geodesic regression. We then use Eq. (3) to give
every hippocampal mesh the same parameterization.

#### Characterization

Consider $q^{\text{start}}$, $q^{\text{end}}$ the meshes corresponding to
hippocampus shapes at the lowest and highest progesterone levels respectively.
Fig. 1 shows that there is a deformation factor of around 10%: each vertex is
displaced by around 10% of the total mesh diameter between the two meshes
shown. Thus, following $\Delta$-Test, we use linear regression to provide a
characterization of 3D shape changes in the hippocampal formation during the
menstrual cycle. Fig. 6 reveals for the first time that the hippocampal
formation shears in response to an increase in progesterone. Additionally,
this computation provides an important educational tool. Clinical
neuroscientists can use the result of our regression model to query, for a
given progesterone level, what is the associated hippocampal shape.

##  6 Conclusion

We have proposed a shape analysis technique to reveal what volumetric analyses
could not: that the overall shape of hippocampal formation changes during
progesterone level fluctuation. The implications for women’s health are
profound. Because each structure of the brain is dedicated to a specific
function, and the hippocampal formation is directly related to functions that
deteriorate in women after menopause, characterizing how the hippocampal
formation changes in response to sex hormones changes is critical. Not only
does it provide a diagnostic for disease prediction, but it also offers a
method to probe relationships between hormone level, hippocampal shape, and
brain health.

Here, we provide a practical method to characterize such changes with slopes
and intercepts learned through approximation schemes for geodesic regression
on the space of 3D discrete surfaces. This work aims to open research avenues
for automated, fast, and statistically sound diagnostics of female brain
health.

## References

  * [1] Martin Bauer, Nicolas Charon, Philipp Harms, and Hsi-Wei Hsieh.  A numerical framework for elastic surface matching, comparison, and interpolation.  International Journal of Computer Vision, 129(8):2425–2444, 2021\.
  * [2] Christopher R Beam, Cody Kaneshiro, Jung Yun Jang, Chandra A Reynolds, Nancy L Pedersen, and Margaret Gatz.  Differences between women and men in incidence rates of dementia and alzheimer’s disease.  Journal of Alzheimer’s disease, 64(4):1077–1083, 2018.
  * [3] M Faisal Beg, Michael I Miller, Alain Trouvé, and Laurent Younes.  Computing large deformation metric mappings via geodesic flows of diffeomorphisms.  International journal of computer vision, 61:139–157, 2005.
  * [4] Janine A Clayton and Francis S Collins.  Policy: Nih to balance sex in cell and animal studies.  Nature, 509(7500):282–283, 2014.
  * [5] AT Du, Nea Schuff, D Amend, MP Laakso, YY Hsu, WJ Jagust, K Yaffe, JH Kramer, B Reed, D Norman, et al.  Magnetic resonance imaging of the entorhinal cortex and hippocampus in mild cognitive impairment and alzheimer’s disease.  Journal of Neurology, Neurosurgery & Psychiatry, 71(4):441–447, 2001.
  * [6] Thomas Fletcher.  Geodesic regression on riemannian manifolds.  In Proceedings of the Third International Workshop on Mathematical Foundations of Computational Anatomy-Geometrical and Statistical Methods for Modelling Biological Shape Variability, pages 75–86, 2011.
  * [7] Liisa AM Galea, Steven R Wainwright, MM Roes, P Duarte-Guterman, C Chow, and DK Hamson.  Sex, hormones and neurogenesis in the hippocampus: hormonal modulation of neurogenesis and potential functional implications.  Journal of neuroendocrinology, 25(11):1039–1061, 2013.
  * [8] Nicolas Guigui, Nina Miolane, Xavier Pennec, et al.  Introduction to riemannian geometry and geometric statistics: from basic theory to implementation with geomstats.  Foundations and Trends® in Machine Learning, 16(3):329–493, 2023.
  * [9] Emmanuel Hartman, Yashil Sukurdeep, Eric Klassen, Nicolas Charon, and Martin Bauer.  Elastic shape analysis of surfaces with second-order sobolev metrics: a comprehensive numerical framework.  International Journal of Computer Vision, 131(5):1183–1209, 2023\.
  * [10] Ian H Jermyn, Sebastian Kurtek, Hamid Laga, Anuj Srivastava, Gerard Medioni, and Sven Dickinson.  Elastic shape analysis of three-dimensional objects.  Springer, 2017.
  * [11] Martin Kilian, Niloy J Mitra, and Helmut Pottmann.  Geometric modeling in shape space.  In ACM SIGGRAPH 2007 papers, pages 64–es. 2007.
  * [12] Sebastian Kurtek, Eric Klassen, Zhaohua Ding, and Anuj Srivastava.  A novel riemannian framework for shape analysis of 3d objects.  In 2010 IEEE computer society conference on computer vision and pattern recognition, pages 1625–1632. IEEE, 2010.
  * [13] Nasim Maleki, Lino Becerra, Jennifer Brawn, Bruce McEwen, Rami Burstein, and David Borsook.  Common hippocampal structural and functional changes in migraine.  Brain Structure and Function, 218:903–912, 2013.
  * [14] Nina Miolane, Nicolas Guigui, Alice Le Brigant, Johan Mathe, Benjamin Hou, Yann Thanwerdas, Stefan Heyder, Olivier Peltre, Niklas Koep, Hadi Zaatiti, Hatem Hajri, Yann Cabanes, Thomas Gerald, Paul Chauchat, Christian Shewmake, Daniel Brooks, Bernhard Kainz, Claire Donnat, Susan Holmes, and Xavier Pennec.  Geomstats: A Python package for Riemannian geometry in machine learning.  Journal of Machine Learning Research, 21:1–9, 2020.
  * [15] Laura Pritschet, Tyler Santander, Caitlin M Taylor, Evan Layher, Shuying Yu, Michael B Miller, Scott T Grafton, and Emily G Jacobs.  Functional reorganization of brain networks across the human menstrual cycle.  Neuroimage, 220:117091, 2020.
  * [16] Lisa R Taxier, Kellie S Gross, and Karyn M Frick.  Oestradiol as a neuromodulator of learning and memory.  Nature Reviews Neuroscience, 21(10):535–550, 2020.
  * [17] Caitlin M Taylor, Laura Pritschet, and Emily G Jacobs.  The scientific body of knowledge–whose body does it serve? a spotlight on oral contraceptives and women’s health factors in neuroimaging.  Frontiers in neuroendocrinology, 60:100874, 2021.
  * [18] Caitlin M Taylor, Laura Pritschet, Rosanna K Olsen, Evan Layher, Tyler Santander, Scott T Grafton, and Emily G Jacobs.  Progesterone shapes medial temporal lobe volume across the human menstrual cycle.  NeuroImage, 220:117125, 2020.
  * [19] P Thomas Fletcher.  Geodesic regression and the theory of least squares on riemannian manifolds.  International journal of computer vision, 105:171–185, 2013.
  * [20] Kjersti Grøtta Vetvik and E Anne MacGregor.  Sex differences in the epidemiology, clinical features, and pathophysiology of migraine.  The Lancet Neurology, 16(1):76–87, 2017.