#  $B_{h}$-sets and Rigidity

Kevin O’Bryant

###### Abstract

A set ${\cal A}$ of nonnegative integers is called a $B_{h}$-set if every
solution to $a_{1}+\dots+a_{h}=b_{1}+\dots+b_{h}$, where $a_{i},b_{i}\in{\cal
A}$, has $\\{a_{1},\dots,a_{h}\\}=\\{b_{1},\dots,b_{h}\\}$ (as multisets). Let
$\gamma_{k}(h)$ be the $k$-th positive element of the greedy $B_{h}$-set.
Previous work has given formulas for $\gamma_{i}(h)$ for $0\leq i\leq 4$. We
give a nontrivial lower bound on $\gamma_{5}(h)$, and a nontrivial upper bound
on $\gamma_{k}(h)$ for $k\geq 5$. Both the upper and lower bounds are of the
form $C_{k}h^{k-1}+O_{k}(h^{k-2})$, with $C_{k}\to 0$ for both.

##  1 Introduction

A set ${\cal A}$ of nonnegative integers is called a $B_{h}$-set if every
solution to

$a_{1}+\dots+a_{h}=b_{1}+\dots+b_{h},\quad a_{i},b_{i}\in{\cal A}$ (1)

has $\\{a_{1},\dots,a_{h}\\}=\\{b_{1},\dots,b_{h}\\}$ (as multisets).
Equivalently, suppose that ${\cal A}=\\{a_{i}:i\in I\\}$, where $I$ may or may
not be infinite. Then ${\cal A}$ is a $B_{h}$-set unless there are distinct
sequences $m_{i},m^{\prime}_{i}$ with $\sum_{i\in I}m_{i}=\sum_{i\in
I}m^{\prime}_{i}=h$ and

$\sum_{i\in I}a_{i}m_{i}=\sum_{i\in I}a_{i}m^{\prime}_{i}.$ (2)

The primary type of literature has been to fix $h$ and to bound the
cardinality $k$ of $B_{h}$-sets in terms of their diameter. For example, if
${\cal A}\subseteq[0,n]$ is a $B_{2h}$-set, then [Green]

$\limsup_{n\to\infty}\frac{|{\cal A}|}{n^{1/(2h)}}\leq\big{(}\sqrt{\pi
h/2}(h/2)!^{2}(1+\epsilon_{h})\big{)}^{1/(2h)},$

with $\epsilon_{h}\to 0$ as $h\to\infty$. Also, there is [BoseChowla] (see
[O'Bryant] for a further development) a $B_{h}$-set ${\cal A}$ with
$\operatorname{diam}({\cal A})<|{\cal A}|^{h}$.

Another area of research concerns the construction and growth of infinite
$B_{h}$-sets. Let ${\cal G}(h)=\\{0=\gamma_{1}<\gamma_{2}<\dots\\}$ be the
lexicographically first infinite $B_{h}$-set, i.e., ${\cal G}(h)$ is the
greedy $B_{h}$-set. For $h=2$, this is called the Mian-Chowla sequence
[MianChowla]; the OEIS [OEIS:mian-chowla] has $25\,000$ terms, the largest of
which is $130\,098\,289\,096$.

The elementary combinatorial bounds on $\gamma_{k}(h)$ are folklore. There are
$\binom{k+h}{k}$ sums with $h$ summands taken from
$\\{\gamma_{0},\dots,\gamma_{k}\\}$, and they are all in the interval
$[0,h\gamma_{k}]$, and so $\gamma_{k}$ cannot be too small. Also,
$\\{\gamma_{0},\dots,\gamma_{k-1},x\\}$ is a $B_{h}$-set for any $x$ except
those of the form

$\frac{1}{r}\left(\gamma_{i_{1}}+\dots+\gamma_{i_{h}}-\gamma_{j_{1}}-\dots-\gamma_{j_{h-r}}\right),\quad
0\leq i,j\leq k-1,1\leq r\leq h.$

Thus, $\gamma_{k}$ is not too large. Quantitatively, these bounds are

$\frac{1}{h}\left[\binom{k+h}{k}-1\right]\leq\gamma_{k}(h)\leq\binom{k+h}{k}^{2}\cdot\frac{kh}{(k+h)^{2}}.$
(3)

For fixed $h$ and $k\to\infty$, line (3) becomes

$\frac{1}{h\cdot h!}\cdot
k^{h}+O_{h}(k^{h-1})\leq\gamma_{k}(h)\leq\frac{h}{h!^{2}}k^{2h-1}+O_{h}(k^{2h-2}).$

There has been no improvement on these bounds.

In this work, we transpose the problem: fix $k$ and bound the $k$-th element
of a $B_{h}$-set as $h$ grows. In this regime, the bounds on line (3) become

$\frac{1}{k!}\cdot
h^{k-1}+O_{k}(h^{k-2})\leq\gamma_{k}(h)\leq\frac{k}{k!^{2}}h^{2k-1}+O_{k}(h^{2k-2}).$

In [Nathanson], the bound $\gamma_{k}(h)\leq h^{k-1}+2h^{k-2}$ is given, which
establishes the order of growth of $\gamma_{k}(h)$ for fixed $k$.

We have the trivial values $\gamma_{0}=0,\gamma_{1}=1,\gamma_{2}=h+1$, and in
[Nathanson], Nathanson shows that $\gamma_{3}=1+h+h^{2}$. In [Nathanson-
OBryant], Nathanson & O’Bryant prove

$\gamma_{4}=\left\lfloor\frac{h+3}{2}\right\rfloor
h^{2}+\left\lfloor\frac{3h+2}{2}\right\rfloor=\frac{1}{2}h^{3}+O(h^{2}).$

The primary result in this work is an improved upper bound $\gamma_{k}$ for
all $k$ (proved in Section 5), and an improved lower bound on $\gamma_{5}$
(proved in Section 4).

######  Theorem 1.1.

Let $\alpha_{5}\coloneqq 0.467214,\alpha_{6}\coloneqq
0.382978,\alpha_{7}\coloneqq 0.269877$, and for $k\geq 7$ set

$\alpha_{k+1}\coloneqq\frac{1}{2^{k}k!}\sum_{j=0}^{k-1}\binom{k-1}{j}\binom{k}{j}2^{j}.$

Then for all $k\geq 5$

$\gamma_{k}\leq\alpha_{k}h^{k-1}+O(h^{k-2}).$

######  Theorem 1.2.

The fifth positive element of the greedy $B_{h}$-set is
$\displaystyle\gamma_{5}\geq\tfrac{1}{8}h^{4}+\tfrac{1}{2}h^{3}$.

#### Maximally Dense ${\mathbf{B}}_{\mathbf{h}}$-sets.

We primarily consider the greedy $B_{h}$-set ${\cal G}(h)$, but also consider
maximally dense $B_{h}$-sets with $k$ elements. Specifically, let $D_{h}(k)$
be the minimum possible value of $a_{k}$ under the constraint that
$0=a_{0}<a_{1}<\dots<a_{k}$ is a $B_{h}$-set. Obviously $D_{h}(0)=0$ and
$D_{h}(1)=1$. We prove in Section 6:

######  Theorem 1.3.

Suppose that $\gcd(a,h+1)=1$. The $3$-element set $\\{0,a_{1},h+1\\}$ is a
$B_{h}$-set of minimum possible diameter; thus $D_{h}(2)=h+1$.

Computations suggest that
$\\{0,1,\frac{1}{2}(h^{2}+h+2),\frac{1}{2}(h^{2}+3h+2)\\}$ is a $B_{h}$-set
with minimum diameter if $h$ is even (if $h$ is odd, increase the last term by
$1$). This suggests that $D_{h}(4)=\frac{1}{2}(h^{2}+3h+2)$ if $h$ is even,
and $D_{h}(4)=\frac{1}{2}(h^{2}+3h+4)$ if $h$ is odd.

#### Rigidity.

This phenomenon, that the transpose problem has polynomial-ish behavior while
the original problem is enigmatic, has recently been noticed in other problems
in additive number theory and is being called _rigidity_. The seminal paper is
[SeniaUlam].

The Ulam set $U(1,n)$ is the set formed by greedily including each positive
integer that is a sum of a pair of smaller terms in exactly one way, and whose
smallest two members are $1$ and $n$. The set $U(1,2)$, the Ulam sequence, is
notoriously difficult to prove anything about, while nevertheless containing
copious apparent qualities. In [SeniaUlam], the following rigidity property is
conjectured: there are integers $m_{i},p_{i},k_{i},r_{i}$ such that for all
$n\geq 4$, the set $U(1,n)\cap[1,Cn]$ is the disjoint union of the intervals
$[m_{i}n+p_{i},k_{i}n+r_{i}]$ (here and throughout this paper, intervals
indicate intervals of integers).

That is, while each $U(1,n)$ is erratic, the transpose problem of describing
$U(1,n)$ as $n\to\infty$ is rigidly polynomial.

##  2 The Greedy ${\mathbf{B}}_{\mathbf{h}}$-set

We set $\gamma_{0}=0$, and inductively set $\gamma_{k+1}$ to be the smallest
integer greater than $\gamma_{k}$ such that
$\\{\gamma_{0},\gamma_{1},\dots,\gamma_{k},\gamma_{k+1}\\}$ is a $B_{h}$-set.
Let ${\cal G}(h)\coloneqq\\{\gamma_{0},\gamma_{1},\dots\\}$ be the _greedy
$B_{h}$-set_.

We define

$\displaystyle F_{r}^{(k)}$
$\displaystyle\coloneqq\left\\{\frac{1}{r}\sum_{i=1}^{k}\gamma_{i}\cdot(m_{i}-m^{\prime}_{i}):m_{i},m^{\prime}_{i}\in\mathbb{N},\sum_{i=1}^{k}m_{i}\leq
h,\sum_{i=1}^{k}m^{\prime}_{i}\leq h-r\right\\},$ (4) $\displaystyle F^{(k)}$
$\displaystyle\coloneqq\bigcup_{r=1}^{h}F_{r}.$ (5)

We draw the reader’s attention to the omission of $\gamma_{0}$ from the
definition of $F_{r}^{(k)}$, and that the sum of multiplicities need not be
$h$. Moreover, in the definition of $F_{r}$ we will often assume (without loss
of generality) that at least one of $m_{i},m^{\prime}_{i}$ is $0$ for each
$i$.

The $\operatorname{mex}$ of a set is the smallest nonnegative integer that is
excluded from the set.

######  Lemma 2.1.

$\displaystyle\gamma_{k+1}=\operatorname{mex}F^{(k)}$.

###### Proof.

If $\gamma_{k}<x<\operatorname{mex}F$, then
$x=\frac{1}{r}\sum_{i=1}^{k}\gamma_{i}\cdot(m_{i}-m^{\prime}_{i})$ as in (4).
Whence, as $\gamma_{0}=0$, we have

$x\cdot
r+\sum_{i=1}^{k}\gamma_{i}m^{\prime}_{i}+\gamma_{0}\cdot(h-r-\sum_{i=1}^{k}m^{\prime}_{i})=\sum_{i=1}^{k}\gamma_{i}m_{i}+\gamma_{0}\cdot(h-\sum_{i=1}^{k}m_{i}),$

proving that $\\{\gamma_{0},\dots,\gamma_{k}\\}\cup\\{x\\}$ is not a
$B_{h}$-set. On the other hand, let $x=\operatorname{mex}F$ and suppose that
$\\{\gamma_{0},\dots,\gamma_{k}\\}\cup\\{x\\}$ is not a $B_{h}$-set, so that
there is a solution to

$xm_{k+1}+\sum_{i=0}^{k}\gamma_{i}m_{i}=xm^{\prime}_{k+1}+\sum_{i=0}^{k}\gamma_{i}m^{\prime}_{i}$

with $\sum_{i=0}^{k+1}m_{i}=\sum_{i=0}^{k+1}m^{\prime}_{i}=h$. As
(inductively) $\\{\gamma_{0},\dots,\gamma_{k}\\}$ is a $B_{h}$-set, it must be
that at least one of $m_{k+1},m^{\prime}_{k+1}$ is positive, say
$m^{\prime}_{k+1}\geq m_{k+1}\geq 0$. If $m^{\prime}_{k+1}=m_{k+1}$, then

$\sum_{i=1}^{k}\gamma_{i}m_{i}+\gamma_{0}\cdot(m_{0}+m_{k+1})=\sum_{i=1}^{k}\gamma_{i}m^{\prime}_{i}+\gamma_{0}\cdot(m^{\prime}_{0}+m^{\prime}_{k+1}),$

contradicting that $\\{\gamma_{0},\dots,\gamma_{k}\\}$ is a $B_{h}$-set. Ergo,
we may set $r\coloneqq m^{\prime}_{k+1}-m_{k+1}>0$. We have

$\sum_{i=1}^{k}\gamma_{i}m_{i}=xr+\sum_{i=1}^{k}\gamma_{i}m^{\prime}_{i},$

which we can rearrange to

$x=\frac{1}{r}\sum_{i=1}^{k}\gamma_{i}\cdot(m_{i}-m^{\prime}_{i}),$

with $\sum_{i=1}^{k}m_{i}=h-m_{0}-m_{k+1}\leq h$ and
$\sum_{i=1}^{k}m^{\prime}_{i}=h-(m^{\prime}_{0}+m^{\prime}_{k+1})\leq h-r$. ∎

##  3 Formulas for ${\mathbf{\gamma}}_{\mathbf{k}}$

The formula $\gamma_{0}=0$ is by definition, and the formula $\gamma_{1}=1$ is
immediate. In [Nathanson], Nathanson provides detailed proofs of the formulas

$\gamma_{2}=h+1,\qquad\text{and}\qquad\gamma_{3}=h^{2}+h+1.$

In [Nathanson-OBryant], Nathanson and the author prove that

$\gamma_{4}=h^{2}\left\lfloor\frac{h+3}{2}\right\rfloor+\left\lfloor\frac{3h+2}{2}\right\rfloor.$

Computation of $\gamma_{5}$ for $h\leq 47$ has yielded a conjectural formula.

######  Conjecture 3.1.

$\gamma_{5}=\frac{h^{4}}{3}+\frac{1}{6}\cdot\begin{cases}5h^{3}+7h^{2}+7h+9,&\text{if
$h\equiv 1\pmod{6}$}\\\ 4h^{3}+10h^{2}+6h+4,&\text{if $h\equiv 2\pmod{6}$}\\\
7h^{3}+5h^{2}+9h+3,&\text{if $h\equiv 3\pmod{6}$}\\\ 406,&\text{if $h=4$}\\\
6h^{3}+6h^{2}+8h+8,&\text{if $h\equiv 4\pmod{6}$ and $h\geq 10$}\\\
6h^{3}+6h^{2}+8h+6,&\text{if $h\equiv 5\pmod{6}$}\\\
5h^{3}+8h^{2}+7h+6,&\text{if $h\equiv 6\pmod{6}$}\\\ \end{cases}$

Paul Voutier [Voutier] has computed $\gamma_{6}$ for $1\leq h\leq 33$,
generating the following sequence:

$h$ | $\gamma_{6}$ | $h$ | $\gamma_{6}$ | $h$ | $\gamma_{6}$
---|---|---|---|---|---
1 | 6 | 13 | 84026 | 25 | 1916949
2 | 30 | 14 | 109870 | 26 | 2361150
3 | 124 | 15 | 156474 | 27 | 2859694
4 | 368 | 16 | 217790 | 28 | 3467661
5 | 926 | 17 | 304910 | 29 | 3989744
6 | 2214 | 18 | 376260 | 30 | 4779270
7 | 4181 | 19 | 510220 | 31 | 5479857
8 | 8043 | 20 | 667130 | 32 | 6449983
9 | 13818 | 21 | 794873 | 33 | 7575912
10 | 23614 | 22 | 1008048 | 34 | ?
11 | 34825 | 23 | 1302947 | 35 | ?
12 | 54011 | 24 | 1629264 | 36 | ?

While all of $\gamma_{0},\dots,\gamma_{4}$ are quasi-polynomials, and
$\gamma_{5}$ appears to be, this author is skeptical that all $\gamma_{k}$
are. However, it is plausible that each $\gamma_{k}$, if $h$ is sufficiently
large, is given by an expression built up from rational functions in $h$ and
floor functions.

Throughout this work, we use the usual interval notation to indicate intervals
of integers; for example $[0,10]\cup[11,20]=[0,20]$, and more generally
$[a,b]\cup[c,d]=[a,d]$ if $a\leq b,c\leq d$ and $a-1\leq c-1\leq b\leq d$.

In an appendix in the arXiv version of this paper, we give the proofs of the
formulas for $\gamma_{2},\gamma_{3},\gamma_{4}$ rewritten in the style of our
proof in the next section that $\gamma_{5}\geq\frac{1}{8}h^{4}-O(h^{3})$.

##  4 A Lower Bound on $\mathbf{\gamma}_{\mathbf{5}}$

We now proceed to prove that
$\gamma_{5}(h)\geq\frac{1}{8}h^{4}+\frac{1}{2}h^{3}$. We assume that $h\geq
5$. The “Iverson Bracket” $\left[\\!\\!\left[\,P\,\right]\\!\\!\right]$ is $1$
if $P$ is true, and is $0$ if $P$ is false. We define functions
$\ell_{i},u_{i}$ for $1\leq i\leq 4$ as follows:

$\begin{array}[]{|rl|rl|}\hline\cr\ell_{4}&:=1&u_{4}&\displaystyle:=\left\lfloor\frac{h+3}{4}\right\rfloor\\\
\ell_{3}(d_{4})&\displaystyle:=d_{4}-\left\lfloor\frac{h+1}{2}\right\rfloor&u_{3}(d_{4})&\displaystyle:=\left\lfloor\frac{h+1}{2}\right\rfloor-
d_{4}\\\
\ell_{2}(d_{4},d_{3})&:=d_{4}+|d_{3}|-h&u_{2}(d_{4},d_{3})&:=h-d_{4}-|d_{3}|\\\
\ell_{1}(d_{4},d_{3},d_{2})&:=\displaystyle
1-h-\sum_{i=1}^{4}d_{i}\left[\\!\\!\left[\,d_{i}\leq
0\,\right]\\!\\!\right]&u_{1}(d_{4},d_{3},d_{2})&\displaystyle:=h-\sum_{i=1}^{4}d_{i}\left[\\!\\!\left[\,d_{i}\geq
0\,\right]\\!\\!\right]\\\ \hline\cr\end{array}$

and now set

$\Delta\coloneqq\bigg{\\{}\left\langle
d_{1},d_{2},d_{3},d_{4}\right\rangle:d_{i}\in\mathbb{Z},\quad\sum_{i=1}^{4}d_{i}\left[\\!\\!\left[\,d_{i}\geq
0\,\right]\\!\\!\right]\leq
h,\quad\sum_{i=1}^{4}d_{i}\left[\\!\\!\left[\,d_{i}<0\,\right]\\!\\!\right]\geq
1-h,\\\ \ell_{4}\leq d_{4}\leq u_{4},\quad\ell_{3}(d_{4})\leq d_{3}\leq
u_{3}(d_{4}),\\\ \ell_{2}(d_{4},d_{3})\leq d_{2}\leq
u_{2}(d_{4},d_{3}),\quad\ell_{1}(d_{4},d_{3},d_{2})\leq d_{1}\leq
u_{1}(d_{4},d_{3},d_{2})\bigg{\\}}.$

The author is unreasonably proud of this definition, in large part because of
the unreasonable amount of time and computation he spent finding it. It is not
dictated by the argument below, and may not be ideal.

For an element $\delta\in\Delta$, its _image_ is the integer
$\delta\cdot\left\langle\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}\right\rangle$.
Clearly the images of elements of $\Delta$ are in $F_{1}$ (defined on line
(4)), and so $\gamma_{5}$ is at least the $\operatorname{mex}$ of the images
of $\Delta$.

We order distinct vectors of the same length by setting $\mu<\nu$ if either
$last(\mu)<last(\nu)$ or both $last(\mu)=last(\nu)$ and
$drop(\mu,-1)<drop(\nu,-1)$. Otherwise, $\mu\geq\nu$. That is, we order
vectors lexicographically considering first their last components. On an
odometer, the units digit climbs from its minimum to its maximum, and then
resets to its minimum at the same moment the tens digit increases by one. If
both the ones and tens digit are at their maximum, then the hundreds digits
increments when the ones and tens digits reset to their minimums. Our ordering
on $\Delta$ is in this spirit, with the least significant digits listed first.

We note that both $\left\langle 0,0,0,1\right\rangle$, with image
$\gamma_{4}$, and $\left\langle
0,h/2,h/4,\lfloor{(h+3)/4}\rfloor\right\rangle$ with image
$\lfloor{(h+3)/4}\rfloor\gamma_{4}\geq\frac{1}{8}h^{4}+\frac{1}{4}h^{3}$, are
in $\Delta$.

Let $\mu_{0}<\mu_{1}<\dots$ be the ordered elements of $\Delta$. We will show
that the images of $\mu_{0},\mu_{1},\dots$, a sequence of integers, increases
by at most $1$ at each term. Thus, the image of $\Delta$ is an interval. As
$\left\langle 0,0,0,1\right\rangle\in\Delta$, we know that $\gamma_{4}$ is an
image, and also

$\left\langle
0,\left\lfloor\frac{h}{2}\right\rfloor,\left\lfloor\frac{h+1}{4}\right\rfloor,\left\lfloor\frac{h+3}{4}\right\rfloor\right\rangle\in\Delta$

which has image

$\left\lfloor\frac{h}{2}\right\rfloor\gamma_{2}+\left\lfloor\frac{h+1}{4}\right\rfloor\gamma_{3}+\left\lfloor\frac{h+3}{4}\right\rfloor\gamma_{4}\geq\frac{1}{8}h^{4}+\frac{1}{2}h^{3}.$

This will prove Theorem 1.2.

Suppose that $\delta<\nu$ are consecutive elements of $\Delta$, say
$\delta=\left\langle d_{1},d_{2},d_{3},d_{4}\right\rangle$ and
$\nu=\left\langle v_{1},v_{2},v_{3},v_{4}\right\rangle$. There are only a four
possibilities:

  1. (_i_)

$d_{1}<u_{1}(d_{4},d_{3},d_{2})$, $\nu=\left\langle
d_{4},d_{3},d_{2},d_{1}+1\right\rangle$;

  2. (_ii_)

$d_{1}=u_{1}(d_{4},d_{3},d_{2}),d_{2}<u_{2}(d_{4},d_{3})$,
$\nu=\left\langle\ell_{1}(d_{4},d_{3},d_{2}+1),d_{2}+1,d_{3},d_{4}\right\rangle$;

  3. (_iii_)

$d_{1}=u_{1}(d_{4},d_{3},d_{2}),d_{2}=u_{2}(d_{4},d_{3}),d_{3}<u_{3}(d_{4})$,
$\nu=\left\langle\ell_{1}(d_{4},d_{3}+1,\ell_{2}(d_{4},d_{3}+1)),\ell_{2}(d_{4},d_{3}+1),d_{3}+1,d_{4}\right\rangle$;

  4. (_iv_)

$d_{1}=u_{1}(d_{4},d_{3},d_{2}),d_{2}=u_{2}(d_{4},d_{3}),d_{3}=u_{3}(d_{4}),d_{4}<u_{4}$,
$\nu=\left\langle\ell_{1}(d_{4}+1,\ell_{3}(d_{4}+1),L),L,\ell_{3}(d_{4}+1),d_{4}+1\right\rangle$,
where $L=\ell_{2}(d_{4}+1,\ell_{3}(d_{4}+1))$.

In each case, the image goes from
$\delta\cdot\left\langle\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}\right\rangle$
to
$\nu\cdot\left\langle\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}\right\rangle$.
Thus, our task is to show that in all cases

$(\nu-\delta)\cdot\left\langle\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}\right\rangle\leq
1.$

Each case is a straightforward usage of the description of the case, the
definitions of $\ell_{i},u_{i}$, and solving algebraic inequalities. One needs
to split case (_ii_) into 4 subcases, depending on the signs of $d_{3},d_{2}$,
and case (_iii_) splits into 2 subcases depending on the sign of $d_{3}$. In
case (_iv_), it is helpful to break into $4$ subcases depending on the residue
of $d_{4}$ modulo $4$.

We show these details in an appendix of the arXiv version of this work.

##  5 An upper bound on ${\mathbf{\gamma}}_{\mathbf{k}}$

In this section we determine a nonincreasing sequence of real numbers
$\alpha_{k}$ with

$\gamma_{k+1}(h)\leq\alpha_{k+1}h^{k}+O(h^{k-1})$

and $\alpha_{k}\to 0$. We will inductively use the value of $\alpha_{k}$ to
give $\alpha_{k+1}$. We take $\alpha_{1}=1$, $\alpha_{2}=1$, $\alpha_{3}=1$,
$\alpha_{4}=\frac{1}{2}$, in accordance with our formulas for $\gamma_{k}$,
$0\leq k\leq 4$, and we assume henceforth that $k\geq 4$ and are working to
bound $\gamma_{k+1}$.

Fix a real number $\beta_{k}$ with
$\frac{1}{2}\alpha_{k}\leq\beta_{k}\leq\alpha_{k}$. We will use

$\beta_{4}\coloneqq 0.406671,\quad\beta_{5}\coloneqq
0.308672,\quad\beta_{6}\coloneqq
0.203975,\quad\beta_{k}\coloneqq\frac{1}{2}\alpha_{k},\quad(k\geq 7).$

The engine of our bound (recall the definitions of $F_{r}^{(k)},F_{r}$ on
lines (4), (5)) is that

$\displaystyle\gamma_{k+1}=\operatorname{mex}F^{(k)}$ $\displaystyle\leq
1+\beta_{k}h^{k}+\left|F^{(k)}\cap(\beta_{k}h^{k},\infty)\right|$
$\displaystyle\leq
1+\beta_{k}h^{k}+\sum_{r=1}^{h}\big{|}F_{r}^{(k)}\cap(\beta_{k}h^{k},\infty)\big{|}.$
(6)

For $r\geq 3$, (and sufficiently large $h$, a hypothesis we use repeatedly in
this section)

$\max F_{r}^{(k)}=\frac{1}{r}\cdot h\gamma_{k}\leq\frac{1}{3}\cdot
h\cdot\left(\alpha_{k}h^{k-1}+O(h^{k-2})\right)<\frac{1}{2}\alpha_{k}h^{k}\leq\beta_{k}h^{k}.$

Thus, for $r\geq 3$ we know that $|F_{r}^{(k)}\cap(\beta_{k}h^{k},\infty)|=0$.
For $r=2$, we arrive at

$\max F_{2}^{(k)}\leq\beta_{k}h^{k}+O(h^{k-1}),$

whence $|F_{2}^{(k)}\cap(\beta_{k}h^{k},\infty)|=O(h^{k-1})$. The bound on
line (6) simplifies to

$\gamma_{k+1}\leq\beta_{k}h^{k}+O(h^{k-1})+|F_{1}^{(k)}\cap(\beta_{k}h^{k},\infty)|.$
(7)

We now stratify $F_{1}^{(k)}$ as $F_{1}(m_{k},m^{\prime}_{k})$, where at least
one of $m_{k},m^{\prime}_{k}$ is $0$, and $0\leq m_{k}\leq h$, and $0\leq
m^{\prime}_{k}\leq h-1$. Specifically

$F_{1}(m_{k},m^{\prime}_{k})\coloneqq(m_{k}-m^{\prime}_{k})\gamma_{k}+\bigg{\\{}\sum_{i=1}^{k-1}(m_{i}-m^{\prime}_{i})\gamma_{i}:\\\
\sum_{i=1}^{k-1}m_{i}\leq h-m_{k},\quad\sum_{i=1}^{k-1}m^{\prime}_{i}\leq
h-1-m^{\prime}_{k},\quad
m_{1}m^{\prime}_{1}=\cdots=m_{k-1}m^{\prime}_{k-1}=0\bigg{\\}}.$

The largest element of $F_{1}(0,m^{\prime}_{k})$ is

$\max F_{1}(0,m^{\prime}_{k})=-m^{\prime}_{k}\gamma_{k}+h\gamma_{k-1}\leq
h\cdot\big{(}\alpha_{k-1}h^{k-2}+O(h^{k-3})\big{)}=\alpha_{k-1}h^{k-1}+O(h^{k-2})<\beta_{k}h^{k},$

(for sufficiently large $h$), so that
$F_{1}(0,m^{\prime}_{k})\cap(\beta_{k}h^{k},\infty)=\emptyset$. We may thus
assume without loss that $m_{k}>0$ and $m^{\prime}_{k}=0$. The largest element
of $F_{1}(m_{k},0)$ is

$\max
F_{1}(m_{k},0)=m_{k}\gamma_{k}+(h-m_{k})\gamma_{k-1}\leq\alpha_{k}\frac{m_{k}}{h}h^{k}+O(h^{k-1}).$

If $m_{k}<h{\beta_{k}}/{\alpha_{k}}$, then $\max F_{1}(m_{k},0)$ is less than
$\beta_{k}h^{k}$ (as always in this section, for large $h$). Ergo, we may
assume that $m_{k}\geq h{\beta_{k}}/{\alpha_{k}}$ (recall our assumption that
$\beta_{k}\leq\alpha_{k}$).

Equation (7) now gives us

$\gamma_{k+1}\leq\beta_{k}h^{k}+O(h^{k-1})+\sum_{\begin{subarray}{c}m_{k}\\\
\frac{\beta_{k}}{\alpha_{k}}h\leq m\leq h\end{subarray}}|F_{1}(m_{k},0)|.$ (8)

Notice that we have dropped the intersection with the interval
$(\beta_{k}h^{k},\infty)$ from our concerns. That might hurt the bound, but
simplicity is its own reward.

We have

$|F_{1}(m_{k},0)|\leq\bigg{|}\bigg{\\{}\big{(}(m_{1},\dots,m_{k-1}),(m^{\prime}_{1},\dots,m^{\prime}_{k-1})\big{)}:\\\
\sum_{i=1}^{k-1}m_{i}\leq h-m_{k},\quad\sum_{i=1}^{k-1}m^{\prime}_{i}\leq
h-1,\quad
m_{1}m^{\prime}_{1}=\cdots=m_{k-1}m^{\prime}_{k-1}=0\bigg{\\}}\bigg{|}.$

Suppose that exactly $j$ of $m^{\prime}_{1},\dots,m^{\prime}_{k-1}$ are
nonzero. By stars-and-bars ($h-1-j$ stars and $j$ bars), there are
$\binom{k-1}{j}\binom{h-1}{h-1-j}$ such tuples. For each such tuple, there are
$h-m_{k}$ stars and $(k-1-j)$ bars in the count of tuples
$(m_{1},\dots,m_{k-1})$; that is, $\binom{h-m_{k}+(k-1-j)}{h-m_{k}}$ valid
tuples $(m_{1},\dots,m_{k-1})$. Altogether, then

$|F_{1}(m_{k},0)|=\sum_{j=0}^{k-1}\binom{k-1}{j}\binom{h-1}{h-1-j}\binom{h-m_{k}+(k-1-j)}{h-m_{k}}.$

Clearly,

$\binom{h-1}{h-1-j}\leq\frac{1}{j!}h^{j}$

and since $h-m_{k}=O(h)$ also

$\displaystyle\binom{h-m_{k}+(k-1-j)}{h-m_{k}}$
$\displaystyle\leq\frac{(h-m_{k}+k-1-j)^{k-1-j}}{(k-1-j)!}$
$\displaystyle=\frac{1}{(k-1-j)!}\sum_{\ell=0}^{k-1-j}\binom{k-1-j}{\ell}(h-m_{k})^{k-1-j-\ell}(k-1-j)^{\ell}$
$\displaystyle=\frac{1}{(k-1-j)!}(h-m_{k})^{k-1-j}+O((h-m_{k})^{k-2-j})$
$\displaystyle=h^{k-1-j}\frac{(1-m_{k}/h)^{k-1-j}}{(k-1-j)!}+O(h^{k-2-j}).$

Thus,

$\displaystyle|F_{1}(m_{k},0)|$
$\displaystyle\leq\sum_{j=0}^{k-1}\binom{k-1}{j}\frac{h^{j}}{j!}\left(\frac{h^{k-1-j}(1-m_{k}/h)^{k-1-j}}{(k-1-j)!}+O(h^{k-2-j})\right)$
$\displaystyle=h^{k-1}\left(\sum_{j=0}^{k-1}\binom{k-1}{j}\frac{(1-m_{k}/h)^{k-1-j}}{j!(k-1-j)!}\right)+O(h^{k-2}).$

Since

$\sum_{\begin{subarray}{c}m_{k}\\\ \frac{\beta_{k}}{\alpha_{k}}h\leq m\leq
h\end{subarray}}(1-m_{k}/h)^{k-1-j}=h\int_{\beta_{k}/\alpha_{k}}^{1}(1-x)^{k-1-j}\,dx+O(1)=h\frac{\left(1-\beta_{k}/\alpha_{k}\right)^{k-j}}{k-j}+O(1),$

we have

$\displaystyle\sum_{\begin{subarray}{c}m_{k}\\\
\frac{\beta_{k}}{\alpha_{k}}h\leq m\leq h\end{subarray}}|F_{1}(m_{k},0)|$
$\displaystyle=h^{k}\left(\sum_{j=0}^{k-1}\binom{k-1}{j}\frac{(1-\beta_{k}/\alpha_{k})^{k-j}}{j!(k-j)!}\right)+O(h^{k-1})$

Line (8) gives us $\gamma_{k+1}\leq\alpha_{k+1}h^{k}+O(h^{k-1})$ provided that

$\beta_{k}+\left(\sum_{j=0}^{k-1}\binom{k-1}{j}\frac{(1-\beta_{k}/\alpha_{k})^{k-j}}{j!(k-j)!}\right)\leq\alpha_{k+1}.$
(9)

Recall our earlier assumption that $\beta_{k}\leq\frac{1}{2}\alpha_{k}$, which
(it turns out) has no slack for $k\geq 7$.

We proceed by choosing $\beta_{k}$ to minimize this expression, giving us the
smallest possible value of $\alpha_{k+1}$.

$\beta_{4}\coloneqq 0.406671,\quad\beta_{5}\coloneqq
0.308672,\quad\beta_{6}\coloneqq
0.203975,\quad\beta_{k}\coloneqq\frac{1}{2}\alpha_{k},\quad(k\geq 7).$

We use $\beta_{k}$ and $\alpha_{k}$ to compute $\alpha_{k+1}$, and for $k\geq
6$ we then use $\alpha_{k+1}$ to compute $\beta_{k+1}$.

$\begin{array}[]{rrl}k&\beta_{k}&\alpha_{k}\\\ 1&&1\\\ 2&&1\\\ 3&&1\\\
4&0.406671&1/2\\\ 5&0.308672&0.467214,\\\ 6&0.203975&0.382978,\\\
7&&0.269877,\\\ \end{array}$

Now set $\beta_{k}=\frac{1}{2}\alpha_{k}$ for $k\geq 7$. We have

$\displaystyle\beta_{k}+\left(\sum_{j=0}^{k-1}\binom{k-1}{j}\frac{(1-\beta_{k}/\alpha_{k})^{k-j}}{j!(k-j)!}\right)$
$\displaystyle=\frac{\alpha_{k}}{2}+\sum_{j=0}^{k-1}\binom{k-1}{j}\frac{2^{-(k-j)}}{j!(k-j)!}$
$\displaystyle=\frac{\alpha_{k}}{2}+\frac{1}{2^{k}k!}\sum_{j=0}^{k-1}\binom{k-1}{j}\binom{k}{j}2^{j}$

Thus, we can set for $k\geq 7$

$\alpha_{k+1}=\frac{\alpha_{k}}{2}+\frac{1}{2^{k}k!}\sum_{j=0}^{k-1}\binom{k-1}{j}\binom{k}{j}2^{j}.$

We comment that

$\alpha_{k+1}\leq\frac{\alpha_{k}}{2}+\frac{1}{2k\cdot
k!}\sum_{j=0}^{k}\binom{k}{j}^{2}=\frac{\alpha_{k}}{2}+\frac{1}{2k\cdot
k!}\binom{2k}{k}\leq\frac{\alpha_{k}}{2}+\frac{4^{k}}{2k\cdot k!},$

and so we have $(2-\epsilon)^{k}\alpha_{k}=o(1)$ for every $\epsilon>0$.

##  6 Optimally Dense ${\mathbf{B}}_{\mathbf{h}}$ sets

We let $D_{h}(k)$ be the minimum possible diameter of a $B_{h}$-set with $k$
elements. Trivially, $D_{h}(1)=0$ and $D_{h}(2)=1$, as $\\{0\\}$ and
$\\{0,1\\}$ are $B_{h}$-sets for every $h$.

Suppose that $0<a<b<h+1$. Since

$ab=\underbrace{a+\dots+a}_{\text{$b$
summands}}+\underbrace{0+\dots+0}_{\text{$h-b$
summands}}=\underbrace{b+\dots+b}_{\text{$a$
summands}}+\underbrace{0+\dots+0}_{\text{$h-a$ summands}}$

the set $\\{0,a,b\\}$ is not a $B_{h}$-set. Thus $F_{h}(3)\geq h+1$.

Suppose that $0<a<h+1$. If $\\{0,a,h+1\\}$ is not a $B_{h}$-set, then there
are nonnegative integers $m_{i},m^{\prime}_{i}$ with
$m_{0}+m_{1}+m_{2}=h,m^{\prime}_{0}+m^{\prime}_{1}=m^{\prime}_{2}=h$ and
$(m_{0},m_{1},m_{2})\not=(m^{\prime}_{1},m^{\prime}_{2},m^{\prime}_{3})$ and

$(h-m_{1}-m_{2})\cdot 0+m_{1}\cdot
a+m_{2}\cdot(h+1)=(h-m^{\prime}_{1}-m^{\prime}_{2})\cdot 0+m^{\prime}_{1}\cdot
a+m^{\prime}_{2}\cdot(h+1).$

This implies that

$(m_{1}-m^{\prime}_{1})a=(m^{\prime}_{2}-m_{2})(h+1),$

giving a common multiple of $a$ and $h+1$. If $a,h+1$ are relatively prime,
then their least common multiple is $a(h+1)$, and so
$|m_{1}-m^{\prime}_{1}|\geq h+1$, which contradicts our supposition that imply
the $m_{i},m^{\prime}_{i}$ are between $0$ and $h$, inclusive. If $a,h+1$ have
a common factor, then $L:=\operatorname{lcm}(a,h+1)<a(h+1)$, whence we may
take $m_{1}=L/a<h+1,m^{\prime}_{1}=0$, $m^{\prime}_{2}=L/(h+1)<a\leq
h,m_{2}=0$, $m_{0}=h-L/a,m^{\prime}_{0}=h-L/(h+1)$.

That is, $\\{0<a<h+1\\}$ is a $B_{h}$-set if and only if $\gcd(a,h+1)=1$.

Experiments indicate that
$\\{0,1,(h^{2}+h+2)/2,(h^{2}+3h+2)/2+\left[\\!\\!\left[\,\text{$h$
odd}\,\right]\\!\\!\right]\\}$ is a $B_{h}$-set with minimum diameter,
suggesting that $D_{h}(4)=(h^{2}+3h+2)/2+\left[\\!\\!\left[\,\text{$h$
odd}\,\right]\\!\\!\right]$.

##  7 Problems That Have Not Been Solved

The first problem is to find and prove formulas for $\gamma_{k}$ for as many
values of $k$ as possible. A useful start to this would be a faster algorithm
for computing ${\cal G}(h)$, or at least a faster implementation.

Lacking formulas for $\gamma_{k}$, we hope for upper bounds superior to that
proved in this work, and for lower bounds applicable for $k\geq 6$ and
superior to that in (3).

All of the known and conjectured formulas for $\gamma_{k}$ are quasi-
polynomials in $h$: there is a modulus $m$ and polynomials $p_{1},\dots,p_{m}$
and $\gamma_{k}(h)=p_{h\bmod m}(h)$. Can one show that all $\gamma_{k}$ are
quasi-polynomials? Moreover, all of the coefficients in all of the quasi-
polynomials are nonnegative. No explanation is known.

Quasi-polynomials arise as Ehrhart polynomials. Is there a more concrete
connection? Is there a region whose Ehrhart polynomial is $\gamma_{4}$ or
$\gamma_{5}$?

In [Nathanson-OBryant], it is noted that we can’t even prove that
$\gamma_{k}(h)<\gamma_{k}(h+1)$.

Of more interest than greedy $B_{h}$-sets are the optimally thick $B_{h}$
sets. Can $D_{h}(k)$, the minimum diameter of a $k$-element $B_{h}$-set, be
computed for $k>3$?

## References

##  Appendix A Formulas for
$\mathbf{\gamma_{2}},\mathbf{\gamma_{3}},\mathbf{\gamma_{4}}$

###  A.1 A Formula for $\mathbf{\gamma_{2}}$

We have $\gamma_{0}=0$ and $\gamma_{1}=1$, and we wish to find $\gamma_{2}$.
By Lemma 2.1,

$\gamma_{2}=\operatorname{mex}\bigg{\\{}\frac{1}{r}(m_{1}-m^{\prime}_{1}):1\leq
r\leq h,\quad m_{1},m^{\prime}_{1}\in\mathbb{N},\\\ m_{1}\leq h,\quad
m^{\prime}_{1}\leq h-r\bigg{\\}}=\operatorname{mex}\\{[-(h-1),h]\\}=h+1.$

###  A.2 A Formula for $\mathbf{\gamma_{3}}$

By Lemma 2.1,

$\gamma_{3}=\operatorname{mex}\bigg{\\{}\frac{1}{r}\big{(}(m_{1}-m^{\prime}_{1})+(h+1)(m_{2}-m^{\prime}_{2})\big{)}:1\leq
r\leq h,\quad m_{1},m^{\prime}_{1},m_{2},m^{\prime}_{2}\in\mathbb{N}\\\
m_{1}+m_{2}\leq h,\quad m^{\prime}_{1}+m^{\prime}_{2}\leq h-r,\quad
m_{1}m^{\prime}_{1}=m_{2}m^{\prime}_{2}=0\bigg{\\}}.$

If $m^{\prime}_{2}>0$, then $m_{2}=0$ and $m_{1}\leq h$, and so

$\frac{1}{r}\big{(}(m_{1}-m^{\prime}_{1})+(h+1)(m_{2}-m^{\prime}_{2})\big{)}=\frac{1}{r}\bigg{(}(m_{1}-m^{\prime}_{1})-(h+1)m^{\prime}_{2}\bigg{)}\leq\frac{1}{r}\bigg{(}h-(h+1)m^{\prime}_{2}\bigg{)}<0.$

Thus we can assume that $m^{\prime}_{2}=0$ and $m_{1}\leq h-m_{2}$. We now
have

$\gamma_{3}=\operatorname{mex}\left\\{\frac{1}{r}\bigg{(}(m_{1}-m^{\prime}_{1})+(h+1)m_{2}\bigg{)}:m_{1}+m_{2}\leq
h,m^{\prime}_{1}\leq h-r,m_{1}m^{\prime}_{1}=0\right\\}.$

With $r=1$, we have $-(h-1)\leq m_{1}-m^{\prime}_{1}\leq h-m_{2}$ and so

$\big{\\{}(m_{1}-m^{\prime}_{1})+(h+1)m_{2}:m_{1}+m_{2}\leq h,\quad
m^{\prime}_{1}\leq h-r,\quad m_{1}m^{\prime}_{1}=0\big{\\}}\\\
=\bigcup_{m_{2}=0}^{h}\bigg{(}(h+1)m_{2}+[-(h-1),h-m_{2}]\bigg{)}.$ (10)

But since the right endpoint of the interval $(h+1)m_{2}+[-(h-1),h-m_{2}]$ is
at least as large as 1 less than the left endpoint of the interval
$(h+1)(m_{2}+1)+[-(h-1),h-(m_{2}+1)]$, the union in (10) is the interval
$[1-h,h(h+1)]$.

If $r\geq 2$, then

$\frac{1}{r}\big{(}m_{1}-m^{\prime}_{1}+(h+1)(m_{2}+m^{\prime}_{2})\big{)}\leq\frac{1}{2}\big{(}h(h+1)\big{)}.$

Thus, $\gamma_{3}=\operatorname{mex}[0,h(h+1)]=h(h+1)+1=h^{2}+h+1.$

###  A.3 A Formula for $\mathbf{\gamma_{4}}$

The formula for $\gamma_{4}$, proved in [Nathanson-OBryant],is more involved,
and the proof requires substantially better organization. The argument given
here is not substantively different from that of [Nathanson-OBryant], although
the exposition is quite different.

The “Iverson Bracket” $\left[\\!\\!\left[\,P\,\right]\\!\\!\right]$ is $1$ if
$P$ is true, and is $0$ if $P$ is false. Set

$\Delta_{r}\coloneqq\left\\{\langle
d_{1},d_{2},d_{3}\rangle:d_{i}\in\mathbb{Z},\quad\sum_{i=1}^{3}d_{i}\left[\\!\\!\left[\,d_{i}>0\,\right]\\!\\!\right]\leq
h,\quad\sum_{i=1}^{3}(-d_{i})\left[\\!\\!\left[\,d_{i}<0\,\right]\\!\\!\right]\leq
h-r\right\\}.$

Then $\gamma_{4}$ is smallest nonnegative integer that is not equal to

$\frac{1}{r}\,\delta\cdot\left\langle\gamma_{1},\gamma_{2},\gamma_{3}\right\rangle$

for any $\delta\in\Delta_{r}$ for any $r\in[1,h]$.

Define

$\Delta\coloneqq\bigg{\\{}\left\langle
m_{1}-m^{\prime}_{1},m_{2}-m^{\prime}_{2},m_{3}\right\rangle:m_{i},m^{\prime}_{i}\in\mathbb{N},\quad\sum_{i=1}^{3}m_{i}\leq
h,\quad\sum_{i=1}^{2}m^{\prime}_{i}\leq h-1,\\\
m_{1}m^{\prime}_{1}=m_{2}m^{\prime}_{2}=m_{3}m^{\prime}_{3}=0,\quad
m^{\prime}_{2}+m_{3}\leq h,\quad m_{3}\leq\frac{h+1}{2}\bigg{\\}},$

which is clearly a subset of $\Delta_{1}$. The cleverness in this proof, if
there is any, is in the choice of the conditions $m^{\prime}_{2}+m_{3}\leq
h,m_{3}\leq\frac{h+1}{2}$, which we don’t explain but are used below.

The values for $m_{3}$ that appear in $\Delta$ are
$[0,\lfloor{(h+1)/2}\rfloor]$. For each value of $m_{3}$, the quantity
$m_{2}-m^{\prime}_{2}$ varies from $-(h-m_{3})$ up to $h-m_{3}$. For each pair
of values $m_{2}-m^{\prime}_{2},m_{3}$, the quantity $m_{1}-m^{\prime}_{1}$
varies through an interval (which depends on whether $m_{2}-m^{\prime}_{2}$ is
nonnegative or negative).

We order the elements of $\Delta$ as follows. We define $\delta=\left\langle
d_{1},d_{2},d_{3}\right\rangle<\delta^{\prime}=\left\langle
d^{\prime}_{1},d^{\prime}_{2},d^{\prime}_{3}\right\rangle$ if there is a $j$
with $d_{3-i}=d^{\prime}_{3-i}$ for $0\leq i<j$ and
$d_{3-j}<d^{\prime}_{3-j}$. In this ordering, we define the _odometer_ of
$\Delta$ to be the numbering $\delta_{0}\coloneqq\left\langle
0,0,0\right\rangle<\delta_{1}<\dots$ of all elements of $\Delta$ that are at
least $\delta_{0}$. Let
$\vec{\gamma}\coloneqq\left\langle\gamma_{1},\gamma_{2},\gamma_{3}\right\rangle=\left\langle
1,h+1,h^{2}+h+1\right\rangle$. The sequence
$\delta_{0}\cdot\vec{\gamma},\delta_{1}\cdot\vec{\gamma},\delta_{2}\cdot\vec{\gamma},\dots$
tends to increase, but need not do so monotonically. We call
$\delta\cdot\vec{\gamma}$ the _image of $\delta$_.

In $\Delta$, the last component $d_{3}$ is never negative. Consecutive
elements of the odometer have either first coordinate going up by $1$, second
coordinate going up by $1$, or third coordinate going up by $1$ (just like a
car odometer, which always has some digit going up by 1, while less
significant digits drop to their minimum possible values). Consecutive
elements of the odometer have one of four shapes:

  1. (_i_)

$\left\langle d_{1},d_{2},d_{3}\right\rangle<\left\langle
d_{1}+1,d_{2},d_{3}\right\rangle$;

  2. (_ii_)

$d_{2}<0$ and $\left\langle d_{1},d_{2},d_{3}\right\rangle=\left\langle
h-m_{3},-m^{\prime}_{2},m_{3}\right\rangle<\left\langle-h+m^{\prime}_{2},-m^{\prime}_{2}+1,m_{3}\right\rangle$;

  3. (_iii_)

$d_{2}\geq 0$ and $\left\langle d_{1},d_{2},d_{3}\right\rangle=\left\langle
h-m_{2}-m_{3},m_{2},m_{3}\right\rangle<\left\langle-(h-1),m_{2}+1,m_{3}\right\rangle$;

  4. (_iv_)

$\left\langle d_{1},d_{2},d_{3}\right\rangle=\left\langle
0,h-m_{3},m_{3}\right\rangle<\left\langle-(h-m_{3}-1),-m_{3},m_{3}+1\right\rangle$.

The restriction $m_{3}\leq(h+1)/2$ in $\Delta$ prevents a fifth shape.

In case (_i_), the image goes up by $1$.

In case (_ii_), the image increases by

$\big{(}\left\langle-h+m^{\prime}_{2},-m^{\prime}_{2}+1,m_{3}\right\rangle-\left\langle
h-m_{3},-m^{\prime}_{2},m_{3}\right\rangle\big{)}\cdot\left\langle
1,h+1,\gamma_{3}\right\rangle=-h+m^{\prime}_{2}+m_{3}+1\leq 1,$

where we have used the condition in the definition of $\Delta$ that
$m^{\prime}_{2}+m_{3}\leq h$. While the “increase” could be negative (or
zero), that is not a concern.

Case (_iii_) can only arise if $\left\langle-(h-1),m_{2}+1,m_{3}\right\rangle$
is in $\Delta$, so that we have $m_{3}+m_{2}+1\leq h$. The image increases by

$\big{(}\left\langle-(h-1),m_{2}+1,m_{3}\right\rangle-\left\langle
h-m_{2}-m_{3},m_{2},m_{3}\right\rangle\big{)}\cdot\left\langle
1,h+1,\gamma_{3}\right\rangle=-h+2+m_{2}+m_{3}\leq 1.$

Case (_iv_) only arises if
$\left\langle-(h-m_{3}-1),-m_{3},m_{3}+1\right\rangle$ is in $\Delta$, so that
we may assume that $m_{3}+1\leq h$. Also, the condition $m_{3}\leq(h+1)/2$ in
$\Delta$ insures that $\left\langle-(h-m_{3}-1),-m_{3},m_{3}+1\right\rangle$
is in $\Delta$. The image increases by

$\big{(}\left\langle-(h-m_{3}-1),-m_{3},m_{3}+1\right\rangle-\left\langle
0,h-m_{3},m_{3}\right\rangle\big{)}\cdot\left\langle
1,h+1,h^{2}+h+1\right\rangle=m_{3}+1-h+1\leq 1.$

Thus, as we move through the odometer, the images start at $0$ and increase by
at most $1$. As $\Delta$ contains $\left\langle
0,\lfloor{\frac{h}{2}}\rfloor,\lfloor{\frac{h+1}{2}}\rfloor\right\rangle$,
which has image
$\lfloor{\frac{h}{2}}\rfloor(h+1)+\lfloor{\frac{h+1}{2}}\rfloor(h^{2}+h+1)$,
we have shown that the image of $\Delta$ contains the interval
$[0,\lfloor{\frac{h}{2}}\rfloor(h+1)+\lfloor{\frac{h+1}{2}}\rfloor(h^{2}+h+1)]$.
And since $\Delta\subseteq\Delta_{1}$, we have shown that the image of
$\Delta_{1}$ contains every natural number in the interval

$\left[0,\left\lfloor\frac{h}{2}\right\rfloor(h+1)+\left\lfloor\frac{h+1}{2}\right\rfloor(h^{2}+h+1)\right].$
(11)

Moreover, $\Delta_{1}$ also contains

$\left\\{\left\langle
d_{1},-\left\lfloor\frac{h+1}{2}\right\rfloor,\left\lfloor\frac{h+3}{2}\right\rfloor\right\rangle:1-h+\left\lfloor\frac{h+1}{2}\right\rfloor\leq
d_{1}\leq h-\left\lfloor\frac{h+3}{2}\right\rfloor\right\\}$

which has image

$\left[\left(h^{2}+h+1\right)\left\lfloor\frac{h+3}{2}\right\rfloor-h\left\lfloor\frac{h+1}{2}\right\rfloor-h+1,h^{2}\left\lfloor\frac{h+3}{2}\right\rfloor+\left\lfloor\frac{3h}{2}\right\rfloor\right].$
(12)

The intervals in (11) and (12) contain every natural number up to and
including $h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor$, which proves
that $\gamma_{4}\geq h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor+1$.

Now suppose that $h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor+1$ is the
image of $\left\langle d_{1},d_{2},d_{3}\right\rangle\in\Delta_{1}$:

$d_{1}+d_{2}(h+1)+d_{3}(h^{2}+h+1)=h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor+1.$
(13)

If $d_{3}\leq(h+1)/2$, then

$d_{1}+d_{2}(h+1)+d_{3}(h^{2}+h+1)\leq(h-\lfloor{\tfrac{h+1}{2}}\rfloor)(h+1)+\lfloor{\tfrac{h+1}{2}}\rfloor(h^{2}+h+1)\\\
<h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor+1.$

If $d_{3}\geq(h+4)/2$, then

$d_{1}+d_{2}(h+1)+d_{3}(h^{2}+h+1)\geq(1-h)(h+1)+\lceil{\tfrac{h+4}{2}}\rceil(h^{2}+h+1)\\\
>h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor+1.$

Thus, $\frac{h+2}{2}\leq d_{3}\leq\frac{h+3}{2}$, i.e.,
$d_{3}=\lfloor{(h+3)/2}\rfloor$. Equation (13) becomes

$d_{2}(h+1)+d_{1}+(h+1)\left\lfloor\frac{h+1}{2}\right\rfloor=\left\lfloor\frac{h}{2}\right\rfloor.$

Reducing this modulo $h+1$ reveals that $d_{1}=\lfloor{h/2}\rfloor$ or
$d_{1}=\lfloor{h/2}\rfloor-(h+1)$, but in the first case
$d_{1}+d_{3}=\lfloor{h/2}\rfloor+\lfloor{(h+3)/2}\rfloor>h$, and so
$\left\langle d_{1},d_{2},d_{3}\right\rangle$ is not in $\Delta_{1}$. In the
second case, we find that $d_{2}=1-\lfloor{(h+1)/2}\rfloor$ and so
$(-d_{1})+(-d_{2})=(h+1)-\lfloor{h/2}\rfloor+\lfloor{(h+1)/2}\rfloor-1>h-1$,
and again $\left\langle d_{1},d_{2},d_{3}\right\rangle\not\in\Delta_{1}$. Thus
no element of $\Delta_{1}$ has
$h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor+1$ has an image.

The largest image of an element of $\Delta_{r}$, for $r\geq 2$, is

$\frac{1}{r}\left\langle
0,0,h\right\rangle\cdot\left\langle\gamma_{1},\gamma_{2},\gamma_{3}\right\rangle\leq\frac{1}{2}h\cdot(h^{2}+h+1)\leq
h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor.$

Thus, $h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor+1$ is not
$\frac{1}{r}\delta\cdot\left\langle\gamma_{1},\gamma_{2},\gamma_{3}\right\rangle$
for any $r$ and any $\delta\in\Delta_{r}$. Thus,
$\gamma_{4}=h^{2}\lfloor{(h+3)/2}\rfloor+\lfloor{3h/2}\rfloor+1.$

##  Appendix B The Details in the Proof that $\gamma_{5}\geq$

We now handle the four possibilities one at a time.

  1. (_i_)

$d_{1}<u_{1}(d_{4},d_{3},d_{2})$, $\nu=\left\langle
d_{4},d_{3},d_{2},d_{1}+1\right\rangle$.

We have
$(\nu-\delta)\cdot\left\langle\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}\right\rangle=\gamma_{1}=1,$
This ends the easiest case.

  1. (_ii_)

$d_{1}=u_{1}(d_{4},d_{3},d_{2}),d_{2}<u_{2}(d_{4},d_{3})$,
$\nu=\left\langle\ell_{1}(d_{4},d_{3},d_{2}+1),d_{2}+1,d_{3},d_{4}\right\rangle$;

We first note that since $\delta\in\Delta$, we know that
$d_{2}\geq\ell_{2}(d_{4},d_{3})$, and since $\nu\in\Delta,$ we know that
$d_{2}+1\leq u_{2}(d_{4},d_{3})$. Combined, we have

$d_{4}+|d_{3}|-h\leq d_{2}\leq h-d_{4}-|d_{3}|-1.$ (14)

We have

$\displaystyle(\nu-\delta)\cdot\left\langle\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}\right\rangle$
$\displaystyle=\big{(}\ell_{1}(d_{4},d_{3},d_{2}+1)-u_{1}(d_{4},d_{3},d_{2})\big{)}\gamma_{1}+\gamma_{2}$
$\displaystyle=\ell_{1}(d_{4},d_{3},d_{2}+1)-u_{1}(d_{4},d_{3},d_{2})+h+1.$

As $d_{4}\geq 1$, we have just $4$ sub-cases to consider based on the sign of
$d_{2}$ and $d_{3}$.

$d_{3}$ | $d_{2}$ | $\ell_{1}(d_{4},d_{3},d_{2}+1)$ | $u_{1}(d_{4},d_{3},d_{2})$
---|---|---|---
negative | negative | $1-h-d_{3}-(d_{2}+1)$ | $h-d_{4}$
negative | nonnegative | $1-h-d_{3}$ | $h-d_{4}-d_{2}$
nonnegative | negative | $1-h-(d_{2}+1)$ | $h-d_{4}-d_{3}$
nonnegative | nonnegative | $1-h$ | $h-d_{4}-d_{3}-d_{2}$
$d_{3}$ | $d_{2}$ | $\ell_{1}(d_{4},d_{3},d_{2}+1)-u_{1}(d_{4},d_{3},d_{2})+h+1$
---|---|---
negative | negative | $1-h+d_{4}-d_{3}-d_{2}$
negative | nonnegative | $2-h+d_{4}-d_{3}+d_{2}$
nonnegative | negative | $1-h+d_{4}+d_{3}-d_{2}$
nonnegative | nonnegative | $2-h+d_{4}+d_{3}+d_{2}$

If $d_{2}<0$, then, we have

$\ell_{1}(d_{4},d_{3},d_{2}+1)-u_{1}(d_{4},d_{3},d_{2})+h+1\leq
1-h+d_{4}+|d_{3}|-d_{2}.$

From line (14), we have $-d_{2}\leq h-d_{4}-|d_{3}|$ and so

$1-h+d_{4}+|d_{3}|-d_{2}\leq 1-h+d_{4}+|d_{3}|+h-d_{4}-|d_{3}|=1.$

Similarly, if $d_{2}\geq 0$ we have

$\ell_{1}(d_{4},d_{3},d_{2}+1)-u_{1}(d_{4},d_{3},d_{2})+h+1\leq
2-h+d_{4}+|d_{3}|+d_{2}\\\ \leq 2-h+d_{4}+|d_{3}|+(h-1-d_{4}-|d_{3}|)=1.$

Thus, in all $4$ sub-cases the image increments by at most $1$, and case
(_ii_) is ended.

  1. (_iii_)

$d_{1}=u_{1}(d_{4},d_{3},d_{2}),d_{2}=u_{2}(d_{4},d_{3}),d_{3}<u_{3}(d_{4})$,
$\nu=\left\langle\ell_{1}(d_{4},d_{3}+1,\ell_{2}(d_{4},d_{3}+1)),\ell_{2}(d_{4},d_{3}+1),d_{3}+1,d_{4}\right\rangle$;

We have

$\displaystyle\delta$ $\displaystyle=\left\langle
u_{1}(d_{4},d_{3},d_{2}),h-d_{4}-|d_{3}|,d_{3},d_{4}\right\rangle$
$\displaystyle\nu$
$\displaystyle=\left\langle\ell_{1}(d_{4},d_{3}+1,d_{4}+|d_{3}+1|-h),d_{4}+|d_{3}+1|-h,d_{3}+1,d_{4}\right\rangle$
$\displaystyle\nu-\delta$
$\displaystyle=\left\langle\ell_{1}(d_{4},d_{3}+1,d_{4}+|d_{3}+1|-h)-u_{1}(d_{4},d_{3},d_{2}),2d_{4}+2\left|d_{3}+\tfrac{1}{2}\right|-2h,1,0\right\rangle$

The image increments by

$\ell_{1}(d_{4},d_{3}+1,d_{4}+|d_{3}+1|-h)-u_{1}(d_{4},d_{3},d_{2})+(2d_{4}+2\left|d_{3}+\tfrac{1}{2}\right|-h)(h+1)+1$
(15)

Moreover, from $\delta\in\Delta$ we know that $d_{3}\geq\ell_{3}(d_{4})$, and
from $\nu\in\Delta$ we know that $d_{3}+1\leq u_{3}(d_{4})$. Thus,

$d_{4}-\left\lfloor\frac{h+1}{2}\right\rfloor\leq
d_{3}\leq\left\lfloor\frac{h+1}{2}\right\rfloor-d_{4}-1.$

We have two subcases depending on whether $d_{3}$ is negative or nonnegative.
Suppose first that $d_{3}<0$. We have

$\displaystyle d_{2}$ $\displaystyle=u_{2}(d_{4},d_{3})=h-d_{4}+d_{3}>0$
$\displaystyle|d_{3}+1|$ $\displaystyle=-d_{3}-1$ $\displaystyle
u_{2}=d_{4}+|d_{3}+1|-h$ $\displaystyle\leq
d_{4}+\left\lfloor\frac{h+1}{2}\right\rfloor-
d_{4}-1-h=-\left\lfloor\frac{h+2}{2}\right\rfloor<0$
$\displaystyle\ell_{1}(d_{4},d_{3}+1,u_{2})$
$\displaystyle=1-h-(d_{3}+1)-(d_{4}-d_{3}-1-h)=1-d_{4}$ $\displaystyle
u_{1}(d_{4},d_{3},d_{2})$
$\displaystyle=h-d_{4}-d_{2}=h-d_{4}-(h-d_{4}+d_{3})=-d_{3}.$

The image increments, substituting into (15), by

$2-d_{4}+d_{3}+(2d_{4}+2\left|d_{3}+\tfrac{1}{2}\right|-h)(h+1).$

Since $d_{3}\leq-1$, we have $2|d_{3}+1/2|=-2d_{3}-1$, giving an increment of

$\displaystyle 2-d_{4}+d_{3}+$ $\displaystyle(2d_{4}-2d_{3}-1-h)(h+1)$
$\displaystyle\leq 2-d_{4}+\left\lfloor\frac{h+1}{2}\right\rfloor-
d_{4}-1+\left(2d_{4}-2d_{4}-2\left\lfloor\frac{h+1}{2}\right\rfloor-1-h\right)(h+1)$
$\displaystyle\leq
1+\left\lfloor\frac{h+1}{2}\right\rfloor-2-\left(2\left\lfloor\frac{h+1}{2}\right\rfloor+1+h\right)(h+1)$
$\displaystyle\leq\frac{h}{2}-(2h+1)(h+1),$

which is negative, for all $h\geq 0$.

Now suppose that $d_{3}\geq 0$.We have

$\displaystyle d_{2}$ $\displaystyle=u_{2}(d_{4},d_{3})=h-d_{4}-d_{3}>0$
$\displaystyle|d_{3}+1|$
$\displaystyle=d_{3}+1\leq\left\lfloor\frac{h+1}{2}\right\rfloor-d_{4}$
$\displaystyle u_{2}=d_{4}+|d_{3}+1|-h$
$\displaystyle\leq\left\lfloor\frac{h+1}{2}\right\rfloor-h<0$
$\displaystyle\ell_{1}(d_{4},d_{3}+1,u_{2})$
$\displaystyle=1-h-u_{2}=-d_{4}-d_{3}$ $\displaystyle
u_{1}(d_{4},d_{3},d_{2})$ $\displaystyle=h-d_{4}-d_{3}-d_{2}.$

The image increments, substituting into (15), by

$\displaystyle 1-d_{4}-d_{3}-$
$\displaystyle(h-d_{4}-d_{3}-d_{2})+(2d_{4}+2d_{3}+1-h)(h+1)+1$
$\displaystyle=2-h+d_{2}+(2d_{4}+2d_{3}+1-h)(h+1)$ $\displaystyle\leq
2-h+h-d_{4}-d_{3}+\bigg{(}2d_{4}+2(\left\lfloor\tfrac{h+1}{2}\right\rfloor-
d_{4}-1)+1-h\bigg{)}(h+1)$
$\displaystyle=2-d_{4}-d_{3}+\bigg{(}2\left\lfloor\frac{h+1}{2}\right\rfloor-1-h\bigg{)}(h+1)$
$\displaystyle\leq 1+0\cdot(h+1)\leq 1.$

Thus, in both sub-cases, the increment is at most $1$, and this ends case
(_iii_).

  1. (_iv_)

$d_{1}=u_{1}(d_{4},d_{3},d_{2}),d_{2}=u_{2}(d_{4},d_{3}),d_{3}=u_{3}(d_{4}),d_{4}<u_{4}$,
$\nu=\left\langle\ell_{1}(d_{4}+1,\ell_{3}(d_{4}+1),L),L,\ell_{3}(d_{4}+1),d_{4}+1\right\rangle$,
where $L:=\ell_{2}(d_{4}+1,\ell_{3}(d_{4}+1))$.

As $d_{4}<u_{4}$, we have $d_{4}\leq\left\lfloor\frac{h-1}{4}\right\rfloor$.
Also, $d_{3}=u_{3}(d_{4})=\left\lfloor\frac{h+1}{2}\right\rfloor-d_{4}$,
$d_{2}=h-\left\lfloor\frac{h+1}{2}\right\rfloor$ and $d_{1}=0$.

Thus,

$\displaystyle\nu$ $\displaystyle=\left\langle-
d_{4},-\left\lfloor\frac{h}{2}\right\rfloor,d_{4}-\left\lfloor\frac{h-1}{2}\right\rfloor,d_{4}+1\right\rangle$
$\displaystyle\delta$ $\displaystyle=\left\langle
0,\left\lfloor\frac{h}{2}\right\rfloor,\left\lfloor\frac{h+1}{2}\right\rfloor-
d_{4},d_{4}\right\rangle$ $\displaystyle\nu-\delta$
$\displaystyle=\left\langle-
d_{4},-2\left\lfloor\frac{h}{2}\right\rfloor,2d_{4}-2\left\lfloor\frac{h+1}{2}\right\rfloor+1,1\right\rangle.$

The increment is then

$-d_{4}-2(h+1)\left\lfloor\frac{h}{2}\right\rfloor+\left(2d_{4}+1-2\left\lfloor\frac{h+1}{2}\right\rfloor\right)(h^{2}+h+1)+\left\lfloor\frac{h+3}{2}\right\rfloor
h^{2}+\left\lfloor\frac{3h+2}{2}\right\rfloor,$

which is more profitably written as

$d_{4}\left(2h^{2}+2h+1\right)-2(h+1)\left\lfloor\frac{h}{2}\right\rfloor+\left(1-2\left\lfloor\frac{h+1}{2}\right\rfloor\right)(h^{2}+h+1)+\left\lfloor\frac{h+3}{2}\right\rfloor
h^{2}+\left\lfloor\frac{3h+2}{2}\right\rfloor.$

We can now use $d_{4}\leq\left\lfloor\frac{h-1}{4}\right\rfloor$. Splitting
into $4$ sub-cases depending the residue of $h$ modulo $4$, we have the
routine problem (for each sub-case) of solving a _quadratic_ inequality.


# A versatile cryogenic system for liquid argon detectors

G. Grauso,11footnotetext: Corresponding author. A. Basco  N. Canci  R. de
Asmundis  F. Di Capua  G. Matteucci  Y. Suvorov  G. Fiorillo

###### Abstract

Detectors for direct dark matter search using noble gases in liquid phase as
detection medium need to be coupled to liquefaction, purification and
recirculation systems. A dedicated cryogenic system has been assembled and
operated at the INFN-Naples cryogenic laboratory with the aim to liquefy and
purify the argon used as active target in liquid argon detectors to study the
scintillation and ionization signals detected by large SiPMs arrays.
The cryogenic system is mainly composed of a double wall cryostat hosting the
detector, a purification stage to reduce the impurities below one part per
billion level, a condenser to liquefy the argon, a recirculation gas panel
connected to the cryostat equipped with a custom gas pump. The main features
of the cryogenic system are reported as well as the performances, long term
operations and stability in terms of the most relevant thermodynamic
parameters.

##  1 Introduction

Noble gases (argon, xenon, neon) in liquid phase are widely used as active
medium in the detectors for direct dark matter search and neutrino
experiments. Argon in particular represents a good sensitive material to this
purpose. In fact liquid argon (LAr) technology is supported at industrial
level to realize large mass detectors; argon provides high scintillation
photon yield222Photon yield $\simeq$40000 $\gamma$/MeV for a minimum ionizing
particle, powerful background rejection through efficient event discrimination
methods with the ionization and scintillation processes produced by particle
interactions, ease of purification and high abundance ($\sim$1% content in the
atmosphere) at reasonable cost [1].
Because of the need for thermodynamic stability and due to the presence of
impurities coming from the detector materials, noble liquid apparatuses need
liquefying, purifying and recirculation cryogenic systems. [2, 3]. Within the
framework of the DarkSide (DS) project a double phase liquid argon Time
Projection Chamber (LAr TPC) detector, hosted in a double wall cryostat, has
been built and installed at INFN-Naples cryogenic laboratory with the aim to
study the scintillation and ionization signals detected by a large array of
SiPMs [4, 5]. To fill, liquify and purify the argon as well as to control and
keep the detector stable in terms of pressure and temperature a dedicated
custom cryogenic system has been built and commissioned.

##  2 Cryogenic system description

The custom cryogenic system is a combined setup made of several subsystems and
can be divided in two segments located respectively inside the clean room and
on the second level of the cryogenic laboratory as shown in fig. 1 [left]. The
clean room hosts the double wall cryostats used for the test setup deployment,
the stainless steel piping gas panel with the valves and the instrumentation
for the argon flow control and monitoring, the gas recirculation pump and the
getter for the argon gas purification. The condenser is placed at the second
level floor above the clean room. This device, which allows to managing two
separated loops for argon and nitrogen, is connected to the cryostat with a
custom-made flexible double-wall transfer line and with a solid cryogenic
transfer line to the external storage liquid nitrogen (LN) tank.

Figure 1: Model of the cryogenic system located on two floors [Left] and
system P&ID [Right].

The working principle of the system is based on two splitted loops for Ar and
N2 as can be inferred by the schematics reported in the Appendix A showing the
P&ID of the system. By opening the manual valve (HW07) and the electric valve
(EV04), the liquid nitrogen coming from the supply tank flows to a small LN
buffer tank inside the condenser where the flow rate of the incoming LN is
regulated by a proportional pneumatic valve (pV10) governed by a proportional
integral derivative (PID) controller set on the nitrogen level. The argon
liquefaction process is made of a pre-cooling phase taking place in two
"brazed plate" heat exchangers, and of the effective liquefaction occurring in
a "tube and shell" heat exchanger. After the heat exchange the nitrogen in gas
phase with a flow rate controlled by a mass flow controller (MFC) is flushed
out into the vent line.

At the beginning of the filling process the gaseous argon (GAr) coming from
the bottle flows to the gas panel, gets directed to the getter and then to the
condenser. The liquified argon drops to the cryostat by gravity. The liquefied
argon in the cryostat evaporates due to the heat exchange with external
ambient and the internal heat load produced by the experimental devices. The
evaporated argon gets pumped through the getter for purification, sent to the
condenser to be liquefied again and reintegrated into the cryostat. The rate
of GAr liquefaction inside the condenser affects the pressure of the Ar in the
cryostat. The cryostat pressure is kept stable by an active PID control
capable of regulating the gaseous nitrogen output flow, adjusting the cooling
power and determining the condenser efficiency.

The condenser has been designed to operate in two different configurations. In
the first configuration, a continuous direct supply of liquid nitrogen is used
as heat exchanger fluid for the argon, while in the second one a cold head
cryocooler can be installed on the dedicated top flange of the condenser to
liquefy the gaseous nitrogen and reintegrate the amount of LN spent during the
argon liquefaction process. The latter configuration would reduce the nitrogen
consumption and would improve the efficiency of the system.

The getter used in the system is the PS$4$-MT$50$-R SAES MonoTorr rare-gas hot
purifier specifically designed to provide ultra-high purity gas for
semiconductor applications [6]. The patented getter alloy operated at elevated
temperatures (up to $300^{\circ}$C) removes impurities by forming irreversible
chemical bonds reducing the outlet impurity levels for O2, H2O, CO, CO2, H2,
N2 and CH4 to < 1 ppb level as required for ionization detection [3, 7, 8, 9].

The cryogenic system can be controlled remotely during operations in a mostly
automatic way. In fact, once the cryostat is filled and the recirculation
started the active components of the cryogenic system are automatically
controlled on the basis of the pressure in the LAr cryostat and of the LN
level in the condenser buffer tank by a Slow Control (SC) system. The SC is
based on a National Instrument LabVIEW Real-Time Application software running
on a CompactRIO controller that collects data from the instrumentation
(pressure trensducers, level transducers, mass flow meter and mass flow
controller) [10].

##  3 Cryogenic system commissioning phase

In the first commissioning test a small 12.9 l double wall vacuum insulated
dewar has been used in order to reduce the volume and filling time, allowing
for a fast feedback on the cryogenic system features. A Metal Bellows gas pump
(Senior Aerospace MB-602) was installed for the argon recirculation. Since the
pump motor speed is fixed (2875 rpm) a manual valve has been placed at the
pump outlet to regulate the flow rate. The cryostat was filled by liquifying
gaseous Ar 6.0 (1 ppm impurities in total) from pressurized bottles. After LAr
filling the recirculation stage was started including the getter in the loop
for argon purification. This commissioning test was mainly focused on the
calibration and adjustment of the several active components of the system in
order to operate with a stable pressure in the cryostat. With this
configuration the minimum pressure fluctuation around the set point of 1100
mbar was found to be lower than about 40 mbar with a standard deviation of 10
mbar, as can be seen in fig. 2.

Figure 2: Cryostat pressure stability trend during the first commissioning
test phase.

##  4 System upgrade and comparative test

After the fine tuning of the control and operational parameters a new test has
been performed with the aim to improve the performance of the system, mainly
in terms of pressure variation after some changes in the system configuration.
For this purpose a larger 300 l cryostat has been used for the test. The
previously used MB 602 gas pump was changed with a tunable speed custom QDrive
recirculation pump. Two expansion tanks were added to reduce the mechanical
vibration and the possibile gas pulsation produced during the recirculation
phase. In addition, to speed up the process the cryostat has been directly
filled with LAr 5.0 purified with an in-line cartridge specifically designed
for liquid argon. To avoid the risk of contamination, the getter was inactive
during the test and excluded from the argon loop. During the whole test period
the system ran at stable conditions obtaining a maximum pressure variation
around the set point lower than 10 mbar, with a standard deviation of 2.4 mbar
as shown in fig. 3. As a comparison, the behaviour of the pressure as a
function of time have been reported in fig. 4 and tab. 1, clearly
demonstrating the reduction of the pressure variation in a range of $\pm$10
mbar around the set point of $\simeq$1100 mbar after the changes in the
cryogenic system.

Figure 3: Behavior of the cryostat pressure stability after the changes in the cryogenic system. Figure 4: Comparison of the pressure trends as function of the time during the performed tests. (Blue: First commissionig test phase, Red: Second commissioning test phase | Average [mbar] | Std. Deviation [mbar] | Max [mbar] | Min [mbar]
---|---|---|---|---
First test phase | 1100 | 10.28 | 1141 | 1060
Second test phase | 1100 | 2.41 | 1110 | 1090
Table 1:  Comparison of the pressure trends during the two test phases

##  5 Conclusions

A versatile custom-made cryogenic system has been installed and commissioned
at INFN Naples cryogenic laboratory aiming to fill in a fast, efficient and
stable way, recirculate and purify the experimental set up with gaseous argon,
and able to keep the cryostat conditions stable over the long period of data-
taking. After the commissioning phase needed for some fine tuning of the
operational parameters and devices, an important improvement of the system was
performed that lead to a sensible reduction of pressure variations in the
cryostat of the order of $\pm$10 mbar over a pressure reference working point
of 1100 mbar. The pressure stability shown during the test allows to keep the
cryogenic parameters of the liquid argon detector stable, and to explore more
accurately the main detector features related to the scintillation and
ionization processes.

##  Appendix A Cryogenic System P&ID

## Acknowledgments

The authors would like to thank Dr. Hanguo Wang for the concept design of the
condenser and the fruitful discussions, Eng. Cary Kendziora, Dr. Alexander
Kish, Prof. Yi Wang, Dr. Xiang Xiao, Dr. Tom Thorpe for the productive
suggestions and indication on the cryogenic system, Eng. Marco Carlini for the
constant help on the engineering design and Gennaro Tortone for the monitoring
system implementation.

## References

  * [1] T. Doke, _Fundamental properties of liquid Argon, Krypton and Xenon as Radiation detector media_ , _Portugal Phys._ 12, 9 (1981).
  * [2] R. Acciarri et al, _Effects of Nitrogen contamination in liquid Argon_ , _JINST_ 5, (2010).
  * [3] R. Acciarri et al, _Oxygen contamination in liquid Argon: combined effects on ionization electron charge and scintillation light_ , _JINST_ 5, (2010).
  * [4] The Global Argon Dark Matter Collaboration, [GADMC], _DarkSide-20k Technical Design Report_ , _DARKSIDE-CSN2-TDR-2112_ (2021).
  * [5] A. Hitachi et al, _Effect of ionization density on time dependence of luminescence from liquid argon and xenon_ , _Phys. Rev. B_ 27, 5279 (1983).
  * [6] SAES Pure Gas, Inc., _PS4-MT50-R Rare Gas Purifier Product Manual_.
  * [7] G. Bakale et al., _Effects of an Electric Field on Electron Attachment to SF 6, N2O and O2 in Liquid Argon and Xenon_, _J. Chem. Phys._ 80 (1976), 2556.
  * [8] Icarus Collaboration, _A study of the electrons drifting over large distances in liquid Argon_ , _Nucl. Inst. and Meth. A_ 275 (1989), 364.
  * [9] E.Buckley et al., _A study of the factors affecting the electron lifetime in ultrapure liquid Argon_ , _Nucl. Inst. and Meth. A_ 305 (1991), 177.
  * [10] www.ni.com.


# Mechanically-driven growth and competition in a Voronoi model of tissues

Louis Brézin<EMAIL_ADDRESS>Kirill S. Korolev<EMAIL_ADDRESS>Department of
Physics, Graduate Program in Bioinformatics and Biological Design Center,
Boston University, Boston, Massachusetts 02215, USA

###### Abstract

The mechanisms leading cells to acquire a fitness advantage and establish
themselves in a population are paramount to understanding the development and
growth of cancer. Although there are many works that study separately either
the evolutionary dynamics or the mechanics of cancer, little has been done to
couple evolutionary dynamics to mechanics. To address this question, we study
a confluent model of tissue using a Self-Propelled Voronoi (SPV) model with
stochastic growth rates that depend on the mechanical variables of the system.
The SPV model is an out-of-equilibrium model of tissue derived from an energy
functional that has a jamming/unjamming transition between solid-like and
liquid-like states. By considering several scenarios of mutants invading a
resident population in both phases, we determine the range of parameters that
confer a fitness advantage and show that the preferred area and perimeter are
the most relevant ones. We find that the liquid-like state is more resistant
to invasion and show that the outcome of the competition can be determined
from the simulation of a non-growing mixture. Moreover, a mean-field
approximation can accurately predict the fate of a mutation affecting
mechanical properties of a cell. Our results can be used to infer evolutionary
dynamics from tissue images, understand cancer-suppressing effects of tissue
mechanics, and even search for mechanics-based therapies.

##  I Introduction

Tissues are communities of cells and undergo evolutionary processes. Cell
competition, in which one cell type takes over another, is particularly
relevant for cancer invasion, but also in developmental processes [1, 2, 3].
Cell competition between a mutant and a resident population is the result of
changes that alter the growth rate of the mutants with respect to the
resident. If a mutant acquires a higher growth rate, it will eventually fixate
and take over the entire population.

In normal tissues, there is a turnover of cells over time [4]: cells undergo
division and are eliminated by apoptosis or cell extrusion [5]. The balance
between cell division and extrusion controls the overall growth of the tissue.
In what follows, we will refer to cell extrusion without distinction between
apoptosis and extrusion because, at our level of description, they both
contribute to the removal of a cell from the monolayer. We are interested in
this paper in the mechanical regulation of growth and its implications in cell
competition [1, 6, 7]. At the cellular level, there are mechanotransductive
proteins that provide mechanical feedback to cells, whose levels can promote
or hinder growth [8, 9, 10]. On a macroscopic level, it has been shown that
exerting pressure on cell spheroids hinders growth, both by lowering the
division rate and by increasing the rate of apoptosis [11]. In cell
monolayers, cells try to maintain a constant homeostatic density by extruding
cells in denser areas, where the pressure is greater [12]. The idea that cell
growth is regulated by pressure is not new [13, 14], but to our knowledge its
implication for cell competition has not been explored outside of continuum
models [15].

To investigate the effect of mechanics on cell competition, we study a
commonly used model of a confluent cell monolayer: the Self-Propelled-Voronoi
(SPV) model. We supplement it with a stochastic growth process with birth and
death rates controlled by tissue mechanics. In the SPV model, the motion of
the cell centers is derived from an energy functional, and the network of
cells is constructed using a Voronoi tessellation. The energy functional is
the same as in the Vertex Model [16, 17], but we model the motion of the cell
centers rather than the vertices. An important aspect of both models is that
they predict a jamming-unjamming transition [18, 19, 20] that has been
suggested to play a central role in several biological processes [21, 22, 23].

After introducing the SPV model for a heterogeneous population composed of two
types of cells and its birth-death dynamics, we simulate invasion scenarios
with all possible single-parameter mutants to determine the outcome of the
competition. We find that the fate of a mutant can be predicted from a
simulation of a mutant-ancestor mixture without cell division or death. Using
this fact and a mean-field approximation, we explain the results of our
simulation by a simple argument based on energy minimization.

##  II Model

Figure 1: Overview of model. a) A Voronoi cell and its parameters. b) Example
of a cell division. The red cell is split in half by a segment that goes
through the center of the cell, with the angle of the segment chosen at
random. Two new cell centers are created at a certain distance normal to this
segment, and a new Voronoi tessellation is computed. c) Example of cell
death/extrusion in which a cell center is removed and a new Voronoi
tessellation is created. d) Initial and e) final configuration of an invasion
simulation, where a small droplet of mutant with a lower preferred area is
introduced with fraction $f=0.1$. A time unit represents 100 simulation time
steps.

###  II.1 Heterogeneous Self-Propelled Voronoi model

We use the Self-Propelled Voronoi (SPV) model to describe the dynamics of a
confluent monolayer of active cells, as found in tissues [19]. In this model,
the cell centers evolve according to equations of motions derived from an
energy functional, with self-propulsion motility. Once the positions of the
cell centers are computed, the shapes of the cells are determined using a
Voronoi tessellation, i.e. by maximizing the size of each cell with respect to
the position of its cell center. Other popular models such as the Vertex Model
[16] or the Cellular Potts model [24] use more degrees of freedom to compute
the shape of cells, but the simple aspect of having the cells themselves as
degrees of freedom make evolutionary dynamics much easier to track, especially
since we have little interest in the details of the cell shapes.

The energy functional for the cells is the same regardless of the chosen
model:

$E=\sum_{\alpha}\frac{K_{\alpha}}{2}(A_{\alpha}-A_{\alpha}^{0})^{2}+\sum_{<i,j>}\Lambda_{ij}l_{ij}+\sum_{\alpha}\frac{\Gamma_{\alpha}}{2}P_{\alpha}^{2},$
(1)

where the index $\alpha$ represents a cell in the tissue. Although the
requirement for the cells to form a Voronoi tessellation appears as an added
constraint compared to the Vertex Model, the dynamics of the network of cells
is very similar in both cases far from the low noise limit [25, 26].

We are interested in the competition between two types of cells, a resident
population labeled with a subscript r and a mutant population labeled with a
subscript m. We assume for now that cells are homogeneous within each type. We
define the total number of cells $N=N_{r}+N_{m}$, and the total energy of the
population of cells is:

$E=\frac{K_{r}}{2}\sum_{\alpha=1}^{N_{r}}(A_{\alpha}-A_{r}^{0})^{2}+\frac{\Gamma_{r}}{2}\sum_{\alpha=1}^{N_{r}}(P_{\alpha}-P_{r}^{0})^{2}+\frac{K_{m}}{2}\sum_{\alpha=N_{r}+1}^{N}(A_{\alpha}-A_{m}^{0})^{2}+\frac{\Gamma_{m}}{2}\sum_{\alpha=N_{r}+1}^{N}(P_{\alpha}-P_{m}^{0})^{2}$
(2)

A common choice is to normalize all lengths by $\sqrt{A_{r}^{0}}$, fixing the
normalized preferred area to unity. However, because our simulation using
cellGPU [27] is setup to have a fixed volume $L^{2}$ set by the initial number
of cells, we believe it is appropriate to normalize the lengths by
$L/\sqrt{N_{0}}$, where $L$ is the length of the simulation box and $N_{0}$ is
the initial number of cells. With this normalization, the average cell area is
determined by the cell number at time $t$:

$\langle
a\rangle=\frac{1}{N(t)}\sum_{\alpha=1}^{N(t)}a_{\alpha}=\frac{N_{0}}{N(t)}$
(3)

The energy is normalized by $K_{r}L^{4}/N_{0}^{2}$. The dimensionless total
energy $e$ reads

$e=\frac{1}{2}\sum_{\alpha=1}^{N_{r}}(a_{\alpha}-a_{r}^{0})^{2}+\frac{\bar{\Gamma}_{r}}{2}\sum_{\alpha=1}^{N_{r}}(p_{\alpha}-p_{r}^{0})^{2}+\frac{\bar{K}}{2}\sum_{\alpha=N_{r}+1}^{N}(a_{\alpha}-a_{m}^{0})^{2}+\frac{\bar{\Gamma}_{m}}{2}\sum_{\alpha=N_{r}+1}^{N}(p_{\alpha}-p_{m}^{0})^{2}$
(4)

There are two dimensionless variables, the areas
$a_{\alpha}=A_{\alpha}N_{0}/L^{2}$ and the perimeters
$p_{\alpha}=P_{\alpha}\sqrt{N_{0}}/L$, and 7 dimensionless parameters. Among
those 7 parameters, four are so-called geometrical parameters
$a_{r}^{0},a_{m}^{0},p_{r}^{0},p_{m}^{0}$ and three are mechanical parameters
$\bar{K},\bar{\Gamma}_{r},\bar{\Gamma}_{m}$.

At each timestep, the force exerted on cell center $i$ is given by

$\bm{F}_{i}=-\frac{\partial
E}{\partial\bm{r}_{i}}=-\sum_{j\in\text{NN}(i)}\frac{\partial
E_{j}}{\partial\bm{r}_{i}}-\frac{\partial E_{i}}{\partial\bm{r}_{i}},$ (5)

where NN$(i)$ refers to the nearest neighbors of cell $i$. The equation of
motion for cell center $i$ is

$\displaystyle\frac{\partial\bm{r}_{i}}{\partial
t}=\mu\bm{F}_{i}+v^{0}\bm{n}_{i}$ (6)

with $\mu$ the mobility of cells taken to be homogeneous amongst cells, and
$\bm{n}_{i}=(\cos\theta_{i},\sin\theta_{i})$ a unit polarization vector that
randomly changes direction over time:

$\partial_{t}\theta_{i}=\eta_{i}(t),\quad\langle\eta_{i}(t)\;\eta_{j}(t^{\prime})\rangle=2D_{r}\delta_{ij}\delta(t-t^{\prime})$
(7)

The cell centers are displaced using a forward Euler discretization of eq. (6)
and the Voronoi tessellation is computed with the new position of the cell
centers.

###  II.2 Coupling between mechanics and growth

While we can infer general trends for the effect of mechanics on growth such
as pressure-induced inhibition of growth [11] or regulation of homeostatic
density [12], pinpointing the exact influence of mechanical variables on
individual cells proves challenging. Experimentally, researchers have access
to macroscopic measures such as the force on a cell assembly or the force
exerted on the substrate using traction force microscopy. However, quantifying
bulk properties such as the stress within a cell monolayer remains elusive.
The stress within a monolayer can be computed using imaging-based inference
methods [28, 29], but since the methods rely on formulating models similar to
the one we are using, it does not provide new information.

In the absence of a clear coupling between mechanics and growth dictated by
experiments, we can think of a generalized linear response of growth to
mechanical changes. If one thinks of the mechanical parameters $a_{0},p_{0}$
in the energy functional as some kind of target homeostatic value for the
tissue, then the birth and death rates of a cell $r_{i}$ can be expressed as

$\displaystyle r^{b}_{i}$
$\displaystyle={r^{b}_{i}}^{0}+\lambda\max\left(0,a-a_{0}+\beta(p-p_{0})\right)$
(8) $\displaystyle r^{d}_{i}$
$\displaystyle={r^{d}_{i}}^{0}-\lambda\min\left(0,a-a_{0}+\beta(p-p_{0})\right)$
(9)

The first term corresponds to the basal components of the birth and death
rates, respectively ${r^{b}_{i}}^{0}$ and ${r^{d}_{i}}^{0}$, which do not
depend on mechanical parameters. The second term is a linear coupling between
the mechanical state of the system and the growth rate $r^{b}_{i}-r^{d}_{i}$,
with a coupling coefficient $\lambda$ and a parameter $\beta$ that weighs the
contribution of the perimeter with respect to the area. The $\min$ and $\max$
in the death and birth rates ensure that each rate is positive. Higher
coupling terms are neglected.

In this first work, we focus on the case $\beta=0$, which means that the
growth rate depends only on the area of the cells, or equivalently on the
density in the tissue. This particular coupling is relevant in epithelial
monolayers, where there is more extrusion of cells in denser regions [12].
This coupling is related to the hydrostatic pressure defined as
$\left.\partial E/\partial A_{i}\right|_{P_{i}}=-K(A_{i}-A_{i}^{0})$. The
hydrostatic pressure represents the normal forces exerted on the edges of the
cells [29].

For simplicity, let us consider the case with the same basal birth and death
rates for all cells, ${r^{b}_{i}}^{0}={r^{d}_{i}}^{0}=r^{0}$. In that case,
the rates are given by:

$\displaystyle r^{b}_{i}$
$\displaystyle=r^{0}+\lambda\max\left(0,a_{i}-a_{0}\right)$ (10)
$\displaystyle r^{d}_{i}$
$\displaystyle=r^{0}-\lambda\min\left(0,a_{i}-a_{0}\right)$ (11)

With our normalization scheme, the average growth rate at time $t$ is fixed by
eq. (3) and depends on the number of cells:

$\langle g_{i}\rangle=\langle r^{b}_{i}-r^{d}_{i}\rangle=\lambda\left(\langle
a_{i}\rangle-a_{0}\right)=\lambda\left(\frac{N(t=0)}{N(t)}-a_{0}\right)$ (12)

The steady state is reached when the average growth rate $\langle
g_{i}\rangle$ vanishes and the steady-state cell number is set by the initial
cell number and the preferred area $N_{SS}=N(t=0)/a_{0}$. We represent the
evolution of cell number over time using different growth couplings in fig. 2.

Another mechanical observable associated with growth is the pressure [11, 30,
14], which corresponds to a special value of $\beta$. The pressure is defined
as the opposite of the isotropic part of the stress tensor, which can be
determined from the forces using the virtual work principle [31]. The
dimensionless pressure $\mathcal{P}$, normalized by $KL^{2}/N_{0}$, is given
by:

$\mathcal{P}=-(a-a^{0})-\frac{\Gamma p}{2a}(p-p^{0}).$ (13)

To linear order in $a-a^{0}$ and $p-p^{0}$, it corresponds to $\beta=\Gamma
p^{0}/(2a^{0})$. The birth and death rates of cell $i$ are expressed as:

$\displaystyle r^{b}_{i}$
$\displaystyle=r^{0}+\lambda\max\left(0,\mathcal{P}_{i}-\mathcal{P}_{i}^{h}\right)$
(14) $\displaystyle r^{d}_{i}$
$\displaystyle=r^{0}-\lambda\min\left(0,\mathcal{P}_{i}-\mathcal{P}_{i}^{h}\right)$
(15)

where $\mathcal{P}^{h}$ is viewed as the homeostatic pressure, since it
determines the value of the pressure in steady state. To keep the number of
initial cells constant, there is a unique value of $\mathcal{P}^{h}$ that
depends on the mechanical parameters.

In a future work, a comparison between different couplings between the growth
rate and mechanical properties and their effect on cell competition will be
studied. Notably, the limit $\beta\to\infty$ is unstable. When $\beta\gg 1$, a
larger $p$ leads to a higher growth rate, regardless of the area. Since there
is no constraint on the perimeter, this leads to infinitely elongated cells.

##  III Results

###  III.1 Growth of an homogeneous tissue

The cell number over time from different couplings between mechanics and
birth-death rates is shown in fig. 2. In all cases,
${r^{b}_{i}}^{0}={r^{d}_{i}}^{0}$ for all cells. We looked at three cases,
$\lambda=0$ (no coupling to mechanics), $\beta=0$ which is density regulation,
and coupling to pressure. In each case, there are five realizations of the
stochastic birth-death process using a Gillespie algorithm. We see that when
there is no coupling to mechanics, the stochastic birth-death process does not
reach a steady state. The cell number is only subject to demographic noise and
there is no effective carrying capacity. However, coupling to mechanical
forces stabilizes the number of cells to a value that depends on the
mechanical parameters. We therefore argue that a mechanical feedback is
necessary to maintain a fixed average cell number when the birth-death process
is noisy.

Figure 2: Left panel: evolution of the number of cells over time for different
regulation of the birth-death process. We notice that purely stochastic
process (blue) leads to large deviations of the cell number over time, while
birth-death processes regulated by mechanics, either through the density (red)
or the pressure (green), lead to a steady-state value of the number of cells.
Right panel: distribution of the time-averaged pressure in the tissue. We
notice that in the solid-like state (solid line) there is a negative pre-
stress value that is not relaxed. However, the liquid-like state (dashed line)
is able to relax this pre-stress.

###  III.2 Invasion by a mutant population

After rescaling physical units, our model contains a total of 15 parameters: 7
energy parameters, 4 growth rate parameters, and 4 motility parameters. Since
we are interested in the competition between mutants and residents, we decided
to fix some parameters of the resident population. We end up with a total of 8
parameters involved in the competition. The parameters and the effect of their
change are recapitulated in table 1.

Parameter Description Effect on fitness of mutant $a_{r,m}^{0}$ Preferred area
Fitness advantage for mutants when $a_{m}^{0}<a_{r}^{0}$. Resident value fixed
to $a_{r}^{0}=1$ $p_{r,m}^{0}$ Preferred perimeter Sets the phase of the
system. Fitness advantage if $p_{r}^{0}>p_{m}^{0}$ $\Gamma_{r,m}$ Perimeter
elasticity No apparent effect. Resident value fixed to $\Gamma_{r}=1$ $K_{m}$
Area elasticity of mutant No apparent effect $\lambda_{r,m}$ Strength of
mechanical coupling to growth No apparent effect. Resident value fixed to
$\lambda_{r}=1$ $r^{0}_{r,m}$ Intrinsic growth rate Provides a fitness
advantage when $r_{m}^{0}>r^{0}_{r}$. Resident value fixed to $r^{0}_{r}=0$
$v_{r,m}^{0}$ Self-advection No effect. Resident value fixed to
$v_{r}^{0}=0.05$ $D_{r,m}$ Rotational diffusion (inverse persistence time) No
effect. Resident value fixed to $D_{r}=1$

Table 1: Simulation parameters and their effect on fitness

The parameters of resident and mutant cells can be thought of as phenotypic
traits that can be changed by mutations. We are interested in the evolutionary
outcome of a change in any of these phenotypic traits. To test this, we
performed numerical invasion experiments in which we introduce a small
fraction $f$ of mutants among a resident population. The mutants have a small
change of the order of $10\%$ in one of the parameters compared to the
residents, and we performed 8 invasion experiments, one for each of the
competition parameters (see table 1). Looking at the evolution of the fraction
$f$ of mutants over time, we can determine if the change is advantageous – $f$
increases on average, eventually leading to fixation – or deleterious – $f$
decreases on average leading to extinction. The evolution of the ratio between
the mutant fraction and the resident fraction is represented in fig. 3 on a
semi-log scale. On average, the mutant fraction evolves according to a
logistic equation $\partial_{t}f=sf(1-f)$, where $s=g_{m}-g_{r}$ is the growth
rate difference between mutants and residents. If $s$ does not depend on the
mutant fraction, the solution to the logistic equation is simply
$\ln(f/(1-f))=st$. Therefore, the straight lines on figure 3 correspond to
growth rate differences that do not depend on the mutant fraction. From figure
3 we identify 3 parameters that provide a fitness difference, while the other
5 parameters do not significantly change the fraction of mutant.
Interestingly, for the 3 parameters that provide a change in fitness, this
change is constant over time, since $\ln(f/(1-f))$ is a linear function of
time.

The effect of the basal rate of the mutant $r_{m}^{0}$ is trivial since it
provides a constant growth rate difference. For the two geometric parameters,
the preferred area $a^{0}$ and the preferred perimeter $p^{0}$, it could be
expected that the growth rate difference generated would change as the mutant
fraction evolves. However, we find that the growth rates of the resident and
mutant population change over time, but the difference between the two remain
constant, mainly due to the volume constraint. Let us take for example the
case where the preferred area is varied, and the mutant has a smaller
preferred area $a_{m}^{0}=0.9$ while the resident preferred area is
$a_{r}^{0}=1$. The volume constraint fixes the average area of the entire
population by the number of cells, according to eq. (3). Initially, the
average area of each cell is fixed to unity. The mutant cells introduced are
trying to be smaller than the resident cells to minimize their energy. With
the average area fixed roughly to unity, the mutant cells will be greater than
their preferred area, but so will the resident cells that grow in size to
compensate for the smaller mutant cells. This effect increases the growth rate
of both mutant and resident cells, depending on their relative fraction and
the total number of cells. However, the difference between the two growth
rates does not depend on the fraction nor the total number of cells. A more
detailed derivation from energetic arguments is provided in section III.3. The
argument for a change in the preferred perimeter $p^{0}$ follows the same
path, since it effectively changes the realized area of the cells. We find
that the effect of the geometric parameters $a^{0}$ and $p^{0}$ on fitness is
stronger in the solid-like phase than it is in the liquid-like phase. We
attribute this effect to the broader distribution of cell areas in the solid-
like phase with respect to the liquid-like phase (see fig. 4).

A small change in the strength of the coupling between the area and the growth
rate $\lambda$ seems to have no effect on the evolution of the population.

Mechanical elasticities $\Gamma,K$ play little role in the fitness when they
are varied alone. However, we show in section III.3 that they control the
magnitude of growth rate difference, even though they cannot create such a
difference.

Figure 3: Evolution of the fraction of mutants from an initial droplet of
fraction $f=0.1$, when the birth and death rates are controlled by the area,
as stated in eqs. 10 and 11. Each panel corresponds to an average of 10
identical simulations, where there is a small change in one parameter for the
mutant population. Solid lines are for simulations where the resident is in a
solid-like state with $p_{r}^{0}=3.3$, and dotted lines are for simulations
where the resident is in a liquid-like state with $p_{r}^{0}=4$

###  III.3 Evolutionary outcome determined by a non-growing mixture

Since the rate of change of the mutant fraction is constant over time, we can
determine the fate of the competition between mutant and resident strains by
simulating a mixture of the two strains at a constant fraction, without
growth. We simulate a 50-50 mixture of mutant and resident strains without
growth and compute the average realized areas for mutant and resident strains,
resp. $\langle a_{m}\rangle$ and $\langle a_{r}\rangle$, which gives the
average growth rate difference between mutant and resident strains
$s=\lambda(\langle a_{m}\rangle-a_{m}^{0}-(\langle a_{r}\rangle-a_{r}^{0}))$.
The results are shown in fig. 5. The sign of the growth rate difference is
well predicted by simulating a mixture without growth, but the actual value of
the growth rate is somewhat different. This discrepancy could be explained by
the fluctuations in cell number induced by the fluctuations in the birth-death
process, since a change in cell number changes the average value of the area
(see eq. 3).

Figure 4: Distribution of areas and perimeter for solid (upper row) and liquid
(lower row) 50-50 mixtures of resident and mutants averaged over time. In the
left column, mutants have a preferred area of $a_{m}^{0}=0.9$ while the one of
residents is $a_{r}^{0}=1$ (horizontal dashed lines). All other parameters are
equal. In the right column, mutants have a preferred perimeter of
$p_{m}^{0}=0.9$ while residents have a preferred perimeter $p_{r}^{0}=1$
(vertical dashed lines). All other parameters are equal and have their
reference value.  Figure 5: Evolution of the fraction of mutants over time
when the preferred area of the mutant is $a_{m}^{0}=0.9$ while the one of
residents is $a_{r}^{0}=1$. The solid line represents the simulation of the
birth-and-death process in a solid state $p_{r}^{0}=p_{m}^{0}=3.63$. All other
parameters are equal and have their reference value. The dotted line is the
estimation from the mean-field approximation given by eq. 19 and the dashed
line is an estimation using the area difference given by a mixture without
growth. We see that although both estimation predict the same evolutionary
outcome, they overestimate the the rate of evolution.

Experimentally, this mean the evolutionary information can be obtained at any
point in time from an heterogeneous population of cells, by averaging the
distribution of areas.

In this section, we present an energy argument to explain why the growth rate
difference between mutants and residents does not depend on the fraction.
Consider a mixture with residents having a preferred area $a_{r}^{0}=1$ and
mutants having a preferred area $a_{m}^{0}$. Let us assume that both residents
and mutants have the same preferred perimeter $p^{0}$. The initial population
is composed of $N_{0}$ cells in a fixed volume $L^{2}$ (an average cell
initially has an area of value 1). We consider that at time $t$ there are
$N_{t}$ cells, a fraction $f=N_{m}/N_{t}$ of mutant cells, and
$1-f=N_{r}/N_{t}$ resident cells. Consider a mean-field approximation to solve
the evolution of identical average cells. In that case, the volume
conservation at time $t$ reads:

$(1-f)a_{r}+fa_{m}=\frac{N_{0}}{N_{t}}$ (16)

Equation (16) sets the average cell size to be $N_{0}/N_{t}$. However, this
value does not minimize the energy of the mixture with two different preferred
areas. We assume that there are small deviations $\delta a_{m},\delta
a_{r}\ll\frac{N_{0}}{N_{t}}$ w.r.t the average value, and write the areas as
$a_{m}=N_{0}/N_{t}+\delta a_{m}$ and $a_{r}=N_{0}/N_{t}+\delta a_{r}$. If the
deviations in areas are small, to first order the perimeters are

$\displaystyle p_{m}\approx p^{0}+c\delta a_{m}$ (17) $\displaystyle
p_{r}\approx p^{0}+c\delta a_{r}$ (18)

We can make the assumption that the proportionality coefficient $c$ is the
same for both mutants and residents since having the same preferred perimeter
means they are similar shaped polygons.

With this perturbation, there is only one degree of freedom, and we can choose
it to be $\delta a_{m}$. By minimizing the energy with respect to $\delta
a_{m}$, we obtain the mean-field equilibrium areas for the mutant and resident
cells. From these areas, we determine the mean-field prediction for the
fitness advantage when the preferred area is varied:

$s_{MF}^{A}=-\lambda\Delta a^{0}\frac{\Gamma c^{2}}{1+\Gamma c^{2}},$ (19)

with $\Delta a^{0}=a_{m}^{0}-a_{r}^{0}$ is the difference in preferred areas.
The sign of $\Delta a^{0}$ determines the evolutionary outcome and supports
the finding from the invasion simulation that $a_{m}^{0}<a_{r}^{0}$ leads to a
fitness advantage.

The exact same thing can be done when the preferred perimeter is varied, and
gives

$s_{MF}^{P}=\lambda\Delta p^{0}\frac{\Gamma c}{1+\Gamma c^{2}},$ (20)

with $\Delta p^{0}=p_{m}^{0}-p_{r}^{0}$ is the difference in preferred
perimeters. Similarly, the sign of $\Delta p^{0}$ determines the evolutionary
outcome and supports the finding from the invasion simulation that
$p_{m}^{0}>p_{r}^{0}$ leads to a fitness advantage.

Both eqs. 19 and 20 confirm that the fitness is independent of frequency. The
growth rate of each type depends on the frequency of the mutant, but the
difference between the growth rate of each does not depend on the frequency.

We can also note that although the normalized perimeter elasticity $\Gamma$
does not play a role without a mutation in the preferred area or perimeter, it
controls the rate of evolution when such a mutation exists. The greater
$\Gamma$ is, the faster evolution occurs. A hand-wavy argument is that a high
perimeter elasticity fixes the shape of the cells, making it harder to deform
in the network to accommodate changes in area.

Both rates depend on the coefficient $c$ that relates a change in area to a
change in perimeter. For a regular n-sided polygon, there is a relationship
$a=1/c_{n}p^{2}$, with $c_{n}=4n\tan(\pi/n)$. In the solid-like phase, when
averaging over time the distribution of areas and perimeters, we indeed
observe a linear relationship between the area and the squared perimeter, as
seen in the first row of fig. 4. Using the coefficient for regular hexagons
$c_{6}\approx 13.9$ fits well the distribution of areas and perimeters in the
solid phase. In the liquid-like phase, the shapes of the cells are not well
described by regular polygons, and there is no a priori relationship between
the area and the perimeter. The simulation of a mixture of resident and mutant
cells shows that there is an area change induced by a change in preferred
perimeter and vice versa. This change in area is smaller than for the solid-
like phase, which explains the smaller rate of change of the mutant fraction
in the liquid-like phase compared to the solid-like phase observed in fig. 3.

##  IV Discussion

Mechanics play a role in the regulation of growth in tissue [12, 11, 30], and
thus can be a driver of cell competition, for example in the case of cancer
invasion [7]. Using an agent-based model of tissue mechanics with
mechanically-regulated growth bridges the gap between evolutionary models of
tumors [32, 33] that are concerned with the genetic diversity in tumors and
macroscopic hydrodynamic or elastic models that look at deformations and
growth induced by heterogeneous cells [34, 35].

Mechanically-regulated growth in a homogeneous tissue controls the average
number of cells in the system even in the presence of demographic noise, by
effectively providing a carrying capacity. To study competition in
heterogeneous tissues, we simulated invasion scenarios by different mechanical
mutants and have identified that the geometrical parameters in the Voronoi
model, the preferred area $a_{0}$ and the preferred perimeter $p_{0}$, confer
a fitness advantage to the mutants, while other parameters do not seem to play
a role. A change in preferred area can be interpreted macroscopically as a
change of homeostatic pressure, where lowering the preferred area is similar
to increasing the homeostatic pressure, i.e. the pressure at which the tissue
is in a homeostatic growth state. In this light, this model agrees with
hydrodynamic models of tumor growth coupled to mechanics [14, 13]. The fitness
advantage, when it exists, is constant over time, which means that it does not
depend on the frequency of the mutant in the tissue. Therefore, we can explain
the long-term evolutionary fate of a mutant – fixation or extinction – either
by energy consideration on a tissue made of average cells or by simulating a
mixture of mutants and residents without growth. This is relevant for
experiments since it means that one does not have to conduct a full
evolutionary experiment to infer the outcome, it is sufficient to know the
average state of tissue at a certain time point.

In this work, we studied in detail the case in which the growth rate depends
on the cell area, and we should note that the results obtained depend on the
chosen mechanical coupling. For example, if one were to consider the birth
rate of a cell to be proportional only to the realized area of the cell
instead of the difference between the realized and preferred areas, a change
in preferred area would have opposite effects on the mutant fraction, as a
larger preferred area would just lead to larger cells. A systematic study in
future works of possible couplings between mechanics and growth could help
discriminate which coupling is most relevant to experimental realizations. We
found that one of the key features of the SPV model, the jamming-unjamming
transition [18, 19, 20], is not a key driver of cell competition. We believe
that there are two reasons why there are no fundamental differences between
the two phases. The first is that the growth process essentially fluidizes the
tissue [36] and gives rise to a non-vanishing effective diffusion, effectively
unjamming the jammed phase. The second reason is that over the time of the
simulation the system is essential well-mixed, there is no sustained spatial
difference between the mutants and the residents.

A next step for this work would be to study the effect of mechanics on genetic
diversity. At larger spatial scales, one would observe nontrivial growth-
induced pressure gradients and traveling fronts, which could change important
aspects of the competition [15, 37, 38].

###### Acknowledgements.

We thank Daniel Sussman for developing cellGPU and making it easily accessible
to the scientific community, and the members of the Theoretical Biophysics
group at BU for helpful discussions and feedback. LB and KK were supported by
NIGMS Grant No. 1R01GM138530-01

## References

  * [1] Y. Tamori and W.-M. Deng, “Cell competition and its implications for development and cancer,” Journal of Genetics and Genomics, vol. 38, pp. 483–495, Oct. 2011.
  * [2] G. Morata and P. Ripoll, “Minutes: Mutants of Drosophila autonomously affecting cell division rate,” Developmental Biology, vol. 42, no. 2, pp. 211–221, 1975.
  * [3] K. S. Korolev, J. B. Xavier, and J. Gore, “Turning ecology and evolution against cancer,” Nature Reviews Cancer, vol. 14, pp. 371–380, May 2014.
  * [4] R. Sender and R. Milo, “The distribution of cellular turnover in the human body,” Nature Medicine, vol. 27, pp. 45–48, Jan. 2021.
  * [5] C. Guillot and T. Lecuit, “Mechanics of Epithelial Tissue Homeostasis and Morphogenesis,” Science, vol. 340, pp. 1185–1189, June 2013.
  * [6] A. Matamoro-Vidal and R. Levayer, “Multiple Influences of Mechanical Forces on Cell Competition,” Current Biology, vol. 29, pp. R762–R774, Aug. 2019.
  * [7] R. Levayer, “Solid stress, competition for space and cancer: The opposing roles of mechanical cell competition in tumour initiation and growth,” Seminars in Cancer Biology, vol. 63, pp. 69–80, June 2020.
  * [8] L. Wagstaff, M. Goschorska, K. Kozyrska, G. Duclos, I. Kucinski, A. Chessel, L. Hampton-O’Neil, C. R. Bradshaw, G. E. Allen, E. L. Rawlins, P. Silberzan, R. E. Carazo Salas, and E. Piddini, “Mechanical cell competition kills cells via induction of lethal p53 levels,” Nature Communications, vol. 7, p. 11373, Apr. 2016.
  * [9] M. Aragona, T. Panciera, A. Manfrin, S. Giulitti, F. Michielin, N. Elvassore, S. Dupont, and S. Piccolo, “A Mechanical Checkpoint Controls Multicellular Growth through YAP/TAZ Regulation by Actin-Processing Factors,” Cell, vol. 154, pp. 1047–1059, Aug. 2013.
  * [10] L. Balasubramaniam, A. Doostmohammadi, T. B. Saw, G. H. N. S. Narayana, R. Mueller, T. Dang, M. Thomas, S. Gupta, S. Sonam, A. S. Yap, Y. Toyama, R.-M. Mège, J. M. Yeomans, and B. Ladoux, “Investigating the nature of active forces in tissues reveals how contractile cells can form extensile monolayers,” Nature Materials, vol. 20, pp. 1156–1166, Aug. 2021.
  * [11] F. Montel, M. Delarue, J. Elgeti, L. Malaquin, M. Basan, T. Risler, B. Cabane, D. Vignjevic, J. Prost, G. Cappello, and J.-F. Joanny, “Stress Clamp Experiments on Multicellular Tumor Spheroids,” Physical Review Letters, vol. 107, p. 188102, Oct. 2011.
  * [12] G. T. Eisenhoffer, P. D. Loftus, M. Yoshigi, H. Otsuna, C.-B. Chien, P. A. Morcos, and J. Rosenblatt, “Crowding induces live cell extrusion to maintain homeostatic cell numbers in epithelia,” Nature, vol. 484, pp. 546–549, Apr. 2012.
  * [13] B. I. Shraiman, “Mechanical feedback as a possible regulator of tissue growth,” Proceedings of the National Academy of Sciences, vol. 102, pp. 3318–3323, Mar. 2005.
  * [14] M. Basan, T. Risler, J. Joanny, X. Sastre‐Garau, and J. Prost, “Homeostatic competition drives tumor growth and metastasis nucleation,” HFSP Journal, vol. 3, pp. 265–272, Aug. 2009.
  * [15] J. Ranft, M. Aliee, J. Prost, F. Jülicher, and J.-F. Joanny, “Mechanically driven interface propagation in biological tissues,” New Journal of Physics, vol. 16, p. 035002, Mar. 2014.
  * [16] R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, “The Influence of Cell Mechanics, Cell-Cell Interactions, and Proliferation on Epithelial Packing,” Current Biology, vol. 17, no. 24, pp. 2095–2104, 2007.
  * [17] S. Alt, P. Ganguly, and G. Salbreux, “Vertex models: from cell mechanics to tissue morphogenesis,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 372, p. 20150520, May 2017.
  * [18] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning, “A density-independent rigidity transition in biological tissues,” Nature Physics, vol. 11, pp. 1074–1079, Dec. 2015.
  * [19] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, “Motility-Driven Glass and Jamming Transitions in Biological Tissues,” Physical Review X, vol. 6, p. 021011, Apr. 2016.
  * [20] M. Chiang and D. Marenduzzo, “Glass transitions in the cellular Potts model,” EPL (Europhysics Letters), vol. 116, p. 28009, Oct. 2016.
  * [21] J. A. Mitchel, A. Das, M. J. O’Sullivan, I. T. Stancil, S. J. DeCamp, S. Koehler, O. H. Ocaña, J. P. Butler, J. J. Fredberg, M. A. Nieto, D. Bi, and J.-A. Park, “In primary airway epithelial cells, the unjamming transition is distinct from the epithelial-to-mesenchymal transition,” Nature Communications, vol. 11, p. 5053, Dec. 2020.
  * [22] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T. Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim, B. Gweon, J. Notbohm, R. Steward Jr, S. Burger, S. H. Randell, A. T. Kho, D. T. Tambe, C. Hardin, S. A. Shore, E. Israel, D. A. Weitz, D. J. Tschumperlin, E. Henske, S. T. Weiss, M. L. Manning, J. P. Butler, J. M. Drazen, and J. J. Fredberg, “Unjamming and cell shape in the asthmatic airway epithelium,” Nature Materials, vol. 14, pp. 1040–1048, Oct. 2015.
  * [23] L. Atia, J. J. Fredberg, N. S. Gov, and A. F. Pegoraro, “Are cell jamming and unjamming essential in tissue development?,” Cells & Development, vol. 168, p. 203727, Dec. 2021.
  * [24] F. Graner and J. A. Glazier, “Simulation of biological cell sorting using a two-dimensional extended Potts model,” Physical Review Letters, vol. 69, pp. 2013–2016, Sept. 1992.
  * [25] D. M. Sussman, M. Paoluzzi, M. Cristina Marchetti, and M. Lisa Manning, “Anomalous glassy dynamics in simple models of dense biological tissue,” EPL (Europhysics Letters), vol. 121, p. 36001, Feb. 2018.
  * [26] D. M. Sussman and M. Merkel, “No unjamming transition in a Voronoi model of biological tissue,” Soft Matter, vol. 14, no. 17, pp. 3397–3403, 2018.
  * [27] D. M. Sussman, “cellGPU: massively parallel simulations of dynamic vertex models,” Computer Physics Communications, vol. 219, pp. 400–406, Oct. 2017.  arXiv: 1702.02939.
  * [28] S. Ishihara and K. Sugimura, “Bayesian inference of force dynamics during morphogenesis,” Journal of Theoretical Biology, vol. 313, pp. 201–211, Nov. 2012.
  * [29] X. Yang, D. Bi, M. Czajkowski, M. Merkel, M. L. Manning, and M. C. Marchetti, “Correlating cell shape and cellular stress in motile confluent tissues,” Proceedings of the National Academy of Sciences, vol. 114, pp. 12663–12668, Nov. 2017.
  * [30] M. Delarue, F. Montel, D. Vignjevic, J. Prost, J.-F. Joanny, and G. Cappello, “Compressive Stress Inhibits Proliferation in Tumor Spheroids through a Volume Limitation,” Biophysical Journal, vol. 107, pp. 1821–1828, Oct. 2014.
  * [31] A. Nestor-Bergmann, G. Goddard, S. Woolner, and O. E. Jensen, “Relating cell shape and mechanical stress in a spatially disordered epithelium using a vertex-based model,” Mathematical Medicine and Biology: A Journal of the IMA, vol. 35, pp. i1–i27, Apr. 2018.
  * [32] I. Bozic, J. G. Reiter, B. Allen, T. Antal, K. Chatterjee, P. Shah, Y. S. Moon, A. Yaqubie, N. Kelly, D. T. Le, E. J. Lipson, P. B. Chapman, L. A. Diaz, B. Vogelstein, and M. A. Nowak, “Evolutionary dynamics of cancer in response to targeted combination therapy,” eLife, vol. 2, p. e00747, June 2013.
  * [33] B. Waclaw, I. Bozic, M. E. Pittman, R. H. Hruban, B. Vogelstein, and M. A. Nowak, “A spatial model predicts that dispersal and cell turnover limit intratumour heterogeneity,” Nature, vol. 525, pp. 261–264, Sept. 2015.
  * [34] A. Taloni, M. Ben Amar, S. Zapperi, and C. A. La Porta, “The role of pressure in cancer growth,” The European Physical Journal Plus, vol. 130, p. 224, Nov. 2015.
  * [35] M. Ben Amar, C. Chatelain, and P. Ciarletta, “Contour Instabilities in Early Tumor Growth Models,” Physical Review Letters, vol. 106, p. 148101, Apr. 2011.
  * [36] J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost, and F. Julicher, “Fluidization of tissues by cell division and apoptosis,” Proceedings of the National Academy of Sciences, vol. 107, pp. 20863–20868, Dec. 2010.
  * [37] K. S. Korolev, M. Avlund, O. Hallatschek, and D. R. Nelson, “Genetic demixing and evolution in linear stepping stone models,” Reviews of Modern Physics, vol. 82, pp. 1691–1718, May 2010.
  * [38] G. Birzu, O. Hallatschek, and K. S. Korolev, “Fluctuations uncover a distinct class of traveling waves,” Proceedings of the National Academy of Sciences, vol. 115, pp. E3645–E3654, Apr. 2018.