# A number-theoretic problem concerning pseudo-real Riemann surfaces

Gareth A. Jones and Alexander K. Zvonkin  School of Mathematical Sciences,
University of Southampton, Southampton SO17 1BJ, UK
<EMAIL_ADDRESS>LaBRI, Université de Bordeaux, 351 Cours de la
Libération, F-33405 Talence Cedex, France<EMAIL_ADDRESS>
###### Abstract.

Motivated by their research on automorphism groups of pseudo-real Riemann
surfaces, Bujalance, Cirre and Conder have conjectured that there are
infinitely many primes $p$ such that $p+2$ has all its prime factors
$q\equiv-1$ mod $(4)$. We use theorems of Landau and Raikov to prove that the
number of integers $n\leq x$ with only such prime factors $q$ is asymptotic to
$cx/\sqrt{\ln x}$ for a specific constant $c=0.4865\ldots$. Heuristic
arguments, following Hardy and Littlewood, then yield a conjecture that the
number of such primes $p\leq x$ is asymptotic to $c^{\prime}\int_{2}^{x}(\ln
t)^{-3/2}dt$ for a constant $c^{\prime}=0.8981\ldots$. The theorem, the
conjecture and a similar conjecture applying the Bateman–Horn Conjecture to
other pseudo-real Riemann surfaces are supported by evidence from extensive
computer searches.

###### Key words and phrases:

Riemann surface, pseudo-real, prime number, Riemann zeta function, Dirichlet
series, Bateman–Horn Conjecture

###### 2010 Mathematics Subject Classification:

11A41, 11M06, 20B25, 30F10

##  1\. Introduction

A compact Riemann surface is real (meaning definable, as a complex algebraic
curve, over $\mathbb{R}$) if and only if it has an orientation-reversing
automorphism of order $2$; it is called pseudo-real if it has an orientation-
reversing automorphism, but none of order $2$. In [2], Bujalance, Cirre and
Conder have proved that for each $g\geq 2$, the maximum order $M^{+}_{\rm
ab}(g)$ of an abelian group of (orientation-preserving) automorphisms of a
pseudo-real Riemann surface of genus $g$ is at least $g$. To demonstrate the
sharpness of this lower bound, their Theorem 4.8 presents a set $\mathcal{A}$
of integers $g$ for which the cyclic group ${\rm C}_{g}$ of order $g$ attains
the upper bound $M^{+}_{\rm ab}(g)=g$. These have the form $g=p+1$ where $p$
is what we will call a BCC prime, defined as follows:

######  Definition 1.1.

A BCC prime is a prime number $p$ such that $p+2$ has only prime factors
$q\equiv-1$ mod $(4)$.

The first sixteen BCC primes are

(1) $5,7,17,19,29,31,41,47,61,67,79,97,101,127,131\;\hbox{and}\;137.$

The authors of [2] show that ${\rm C}_{g}$ is the orientation-preserving
subgroup of ${\rm C}_{2g}$ in its action on a compact Riemann surface
$\mathcal{S}$ of genus $g$, with ${\mathcal{S}}/{\rm C}_{2g}$ having signature
$(1;-;[2,2,g];\\{-\\})$. (This is both a description of the quotient space as
an orbifold, namely the real projective plane with cone points of orders $2,2$
and $g$, and an encoded presentation of the non-euclidean crystallographic
group which uniformises it.) If $g\in{\mathcal{A}}$ then $g$ is even, so the
only involution in ${\rm Aut}\,{\mathcal{S}}$ preserves orientation, and hence
$\mathcal{S}$ is pseudo-real. They then use a detailed case-by-case argument
to show that if also $g>30$ then no pseudo-real Riemann surface of genus $g$
has a larger abelian group of automorphisms. For example, since $p+2$ has no
prime factors $q\equiv 1$ mod $(4)$, ${\rm Aut}\,{\rm C}_{p+2}$ (of order
$\phi(p+2)$) has no elements of order $4$, a fact used in subcase (2c) of the
proof of Theorem 4.8 to eliminate ${\rm C}_{g+1}={\rm C}_{p+2}$. Although our
main emphasis here is on number theory and computation, for readers interested
in Riemann surfaces we have described the construction of the surfaces
$\mathcal{S}$ in more detail in an appendix, Section 10.

In [2], $\mathcal{A}$ is described as ‘a very large and possibly infinite
set’. If $\mathcal{A}$ is finite, then conceivably there is a better lower
bound for $M^{+}_{\rm ab}(g)$, valid for all sufficiently large $g$ (in
particular, larger than all those in $\mathcal{A}$, which can then be regarded
as small exceptions). To avoid this possibility, it is important to know
whether $\mathcal{A}$ is infinite. We are therefore interested in the
following problem:

######  Problem 1.2.

Are there infinitely many BCC primes, or equivalently, is the set
$\mathcal{A}$ infinite?

Of course, the answer will be positive if one can prove that there are
infinitely many prime pairs $p,p+2$ with $p\equiv 1$ mod $(4)$, but even
without this extra congruence condition, the existence or otherwise of
infinitely many prime pairs is a very difficult open problem. Nevertheless,
this would follow from a proof of various conjectures in Number Theory, such
as the Bateman–Horn Conjecture [1] or Schinzel’s Hypothesis H [15]. Here we
give a heuristic argument, supported by computational evidence for their
asymptotic distribution, for the following

######  Conjecture 1.3.

There are infinitely many BCC primes.

(It follows from a theorem of Iwaniec [7] that there are infinitely many
primes $p$ such that $p+2$ has only prime factors $q\equiv 1$ mod $(4)$,
whereas the corresponding result with $q\equiv-1$ mod $(4)$ is unproved; see
the Digression in Section 3 for details.)

In support of Conjecture 1.3 we use theorems of Landau [8] and of Raikov [13]
to give, in Theorem 3.1, an asymptotic estimate of the form $cx/\sqrt{\ln x}$
(with $c$ a specific constant) for the number of integers $n\leq x$ with only
prime factors $q\equiv-1$ mod $(4)$. Heuristic arguments, inspired by Hardy
and Littlewood [6], then allow us to make Conjecture 6.1, that the number of
BCC primes $p\leq x$ is asymptotic to

$c^{\prime}\negthickspace\int_{2}^{x}\negthickspace\frac{dt}{(\ln t)^{3/2}}$

for a specific constant $c^{\prime}$. We give numerical evidence, based on
computer searches, to support both Theorem 3.1 and Conjecture 6.1. In
addition, in Section 9, we briefly consider another set of primes, defined by
similar but more complicated conditions, also arising from a construction in
[2].

The formulae we obtain involve various multiplicative factors. These are
defined in terms of well-known constants such as $e$, $\gamma$ and $\pi$,
together with constants such as $C(k,u)$ (defined later) which have been
computed elsewhere to over $100$ decimal places, so in principle they can be
computed with similar accuracy. However, for simplicity we have in most cases
presented numerical data to about ten significant figures, since this is
adequate for the arguments we wish to present.

##  2\. Raikov’s Theorem

If $\mathcal{P}$ is a set of prime numbers, then let us define an integer
$n\in\mathbb{N}$ to be a $\mathcal{P}$-integer if all its prime factors are
elements of $\mathcal{P}$, and let ${\mathcal{P}}^{*}$ denote the set of all
$\mathcal{P}$-integers. If we define $g_{n}=1$ or $0$ as
$n\in{\mathcal{P}}^{*}$ or not, then the function

(2) $g(x):=\sum_{n\leq x}g_{n}$

gives the number of $\mathcal{P}$-integers $n\leq x$. We will consider the
asymptotic behaviour of $g(x)$ as $x\to\infty$.

In [8] Landau showed that if, for some $k\in{\mathbb{N}}$, $\mathcal{P}$ is
the set of all primes in the union of $l$ distinct congruence classes of units
mod $(k)$, then

$g(x)\sim\frac{cx}{(\ln x)^{1-l/\phi(k)}}\quad\hbox{as}\quad x\to\infty$

for some constant $c>0$. In the context of Problem 1.2, if we take $k=4$ and
define ${\mathcal{P}}$ to be the set of primes in the congruence class
$[-1]\in{\mathbb{Z}}^{*}_{4}$, we see that

$g(x)\sim\frac{cx}{\sqrt{\ln x}}\quad\hbox{as}\quad x\to\infty$

for some $c>0$. Our aim is to determine this constant $c$, and then to compare
the resulting estimates with the actual values of $g(x)$ for various $x$.

For any set $\mathcal{P}$ of primes, the function $n\mapsto g_{n}$ is
completely multiplicative, so the corresponding Dirichlet series

(3) $F(s):=\sum_{n=1}^{\infty}\frac{g_{n}}{n^{s}}$

has an Euler product expansion

(4) $F(s)=\prod_{q\in{\mathcal{P}}}\left(1-\frac{1}{q^{s}}\right)^{-1}.$

By comparison with the Riemann zeta function

(5) $\zeta(s):=\sum_{n=1}^{\infty}\frac{1}{n^{s}}=\prod_{q\;{\rm
prime}}\left(1-\frac{1}{q^{s}}\right)^{-1},$

this series and product for $F(s)$ converge for all $s\in\mathbb{C}$ with
${\rm Re}(s)>1$.

A theorem of Raikov [13], as given in [3, Theorem 2.4.1], states the
following:

######  Theorem 2.1.

Let $F(s)=\sum_{n\geq 1}a_{n}/n^{s}$ be a Dirichlet series with non-negative
coefficients, converging for ${\rm Re}(s)>1$. Suppose that $F(s)$ extends
analytically at all points on ${\rm Re}(s)=1$ apart from $s=1$, and that at
$s=1$ we can write

(6) $F(s)=\frac{H(s)}{(s-1)^{1-\alpha}}$

for some $\alpha\in{\mathbb{R}}$ and some $H(s)$ holomorphic in the region
${\rm Re}(s)\geq 1$ and nonzero there. Then

(7) $\sum_{n\leq x}a_{n}\sim\frac{cx}{(\ln x)^{\alpha}}$

as $x\to\infty$, with

(8) $c=\frac{H(1)}{\Gamma(1-\alpha)}$

where $\Gamma$ is the Gamma function

$\Gamma(z)=\int_{0}^{\infty}\negthickspace t^{z-1}e^{-t}\,dt$

for ${\rm Re}(z)>0$.

In simple terms, (6) means that $F(s)$ has a pole of order $1-\alpha$ at
$s=1$, in the sense that it is a non-zero holomorphic multiple of
$1/(s-1)^{1-\alpha}$ near $s=1$, and $H(1)$ is the limiting value of the
multiplier as $s\to 1$.

In the next section we will determine the function $H$ and hence the constants
$\alpha$ and $c$ for the sequence $(g_{n})$. We will do this by expressing the
corresponding function $F$ in terms of $\zeta$ and related functions, so that
the required analyticity conditions for $F$ will follow from those of $\zeta$.

##  3\. Distribution of $\mathcal{P}$-integers

In the case of Problem 1.2 we take

${\mathcal{P}}:=\\{q\mid q\;\hbox{is prime and}\;q\equiv-1\;{\rm
mod}\,(4)\\}.$

Thus the $\mathcal{P}$-integers $n\leq 100$ are

$1,3,7,9,11,19,21,23,27,31,33,43,47,49,57,59,63,67,69,71,77,79,81,83,93,99.$

Note that we include $1$ as the empty product of primes $q\equiv-1$ mod $(4)$.
These $26$ integers represent just over half of the odd integers $n\leq 100$.
However, this proportion decreases towards $0$ as the upper bound increases:
see Theorem 3.1 and Table 1, later in this section.

In order to determine the corresponding constants $c$ and $\alpha$ we can
determine the right-hand side of (6) by comparing $F(s)$, or more precisely
$F(s)^{2}$, with $\zeta(s)$, which is meromorphic with a simple pole of
residue $1$ at $s=1$ (see [4, Theorem 8.28], for example).

For ${\rm Re}(s)>1$ we have

$\frac{F(s)^{2}}{\zeta(s)}=\left(1-\frac{1}{2^{s}}\right)\cdot\prod_{q\equiv-1(4)}\left(1-\frac{1}{q^{s}}\right)^{-2}\cdot\prod_{q\equiv\pm
1(4)}\left(1-\frac{1}{q^{s}}\right)=\left(1-\frac{1}{2^{s}}\right)\cdot\prod_{q\equiv\pm
1(4)}\left(1-\frac{1}{q^{s}}\right)^{\pm 1}.$

so that

(9) $(s-1)F(s)^{2}=\left(1-\frac{1}{2^{s}}\right)\cdot\prod_{q\equiv\pm
1(4)}\left(1-\frac{1}{q^{s}}\right)^{\pm 1}\cdot(s-1)\zeta(s).$

As $s\to 1$ we have

$\left(1-\frac{1}{2^{s}}\right)\to\frac{1}{2},\quad\prod_{q\equiv\pm
1(4)}\left(1-\frac{1}{q^{s}}\right)^{\pm 1}\to\prod_{q\equiv\pm
1(4)}\left(1-\frac{1}{q}\right)^{\pm 1}\quad\hbox{and}\quad(s-1)\zeta(s)\to
1.$

In order to deal with the second of these limits, we can use some results of
Languasco and Zaccagnini [9, 10, 11]. For each integer $k\geq 3$ and each unit
$u$ mod $(k)$ they define a non-zero Mertens-type constant $C(k,u)$ by the
asymptotic estimate

(10) $P(k,u;x):=\prod_{x\geq q\equiv
u(k)}\left(1-\frac{1}{q}\right)=\frac{C(k,u)}{(\ln
x)^{1/\phi(k)}}+O\left(\frac{1}{(\ln
x)^{1+1/\phi(k)}}\right)\quad\hbox{as}\quad x\to\infty,$

where the product is over all primes $q\equiv u$ mod $(k)$ such that $q\leq
x$, and $\phi$ is Euler’s totient function. In [10, Equation (2)] they show
that

(11) $C(k,u)^{\phi(k)}=e^{-\gamma}\prod_{q\,\rm
prime}\left(1-\frac{1}{q}\right)^{\alpha(q)}$

where $\gamma$ is the Euler–Mascheroni constant, the product is now over all
primes $q$, and $\alpha(q):=\phi(k)-1$ or $-1$ as $q\equiv u$ mod $(k)$ or
not. By taking $k=4$ and $u=1$ in (11) we see that

$\prod_{q\equiv\pm 1(4)}\left(1-\frac{1}{q}\right)^{\pm
1}=e^{\gamma}C^{2}\left(1-\frac{1}{2}\right)=\frac{e^{\gamma}C^{2}}{2}$

where $C:=C(4,1)$.

It therefore follows by taking limits as $s\to 1$ in (9) that

$H(1)^{2}=\lim_{s\to 1}\,(s-1)F(s)^{2}=\frac{e^{\gamma}}{4}C^{2}\neq 0,$

so that this particular instance of (6) becomes $F(s)=H(s)/(s-1)^{1/2}$ and
hence

$\alpha=\frac{1}{2}\quad\hbox{and}\quad
H(1)=\sqrt{e^{\gamma}}\cdot\frac{C}{2}.$

Since $\Gamma(\frac{1}{2})=\sqrt{\pi}$ we therefore have

(12)
$c=\frac{H(1)}{\Gamma(\frac{1}{2})}=\sqrt{\frac{e^{\gamma}}{\pi}}\cdot\frac{C}{2},$

giving the following result111Here we cannot resist expressing our pleasure at
the appearance of the constant $\sqrt{e^{\gamma}/\pi}$, in which the three
basic constants $e$, $\pi$ and $\gamma$ of Analysis are united by three of the
basic operations (division, exponentiation and taking square roots) of
Algebra; its application to Number Theory adds to the pleasure.:

######  Theorem 3.1.

The function $g(x)$ satisfies

(13) $g(x)\sim\frac{cx}{\sqrt{\ln
x}}\;=\sqrt{\frac{e^{\gamma}}{\pi}}\cdot\frac{C}{2}\cdot\frac{x}{\sqrt{\ln
x}}\quad\hbox{as}\quad x\to\infty.$

In [10, 11] Languasco and Zaccagnini have evaluated many of these constants
$C(k,u)$ to over $100$ decimal places. These include

(14) $C=C(4,1)=1.2923041571286886071091383898704320653429\ldots,$

which, together with

$\sqrt{\frac{e^{\gamma}}{\pi}}=0.7529495060464205959354997575876048108386\ldots,$

allows us to evaluate

(15) $c=0.4865198883858909971272456405868234055382\ldots$

We use this in Table 1 to test the accuracy of the resulting estimates given
by Theorem 3.1. The second column gives the number $g(x)$ of
${\mathcal{P}}$-integers $n\leq x$ for $x=10^{k}$ with $k=1,\ldots,9$. For
example, the $25$ such integers $n\leq 10^{2}$ are

$3,\,7,\,9,\,11,\,19,\,21,\,23,\,27,\,31,\,33,\,43,\,47,\,49,\,57,\,59,\,63,\,67,\,69,\,71,\,77,\,79,\,81,\,83,\,93,\,99.$

Just $12$ of these correspond to BCC primes $p=n-2$, namely those primes
displayed in Section 1, with $12$ associated genera
$g=p+1=n-1\in{\mathcal{A}}$. The third column of Table 1 gives the estimates
$cx/\sqrt{\ln x}$ for $g(x)$, rounded to the nearest integer, and the fourth
column gives their errors.

$x$ | $g(x)$ | $cx/\sqrt{\ln x}$ | error | $E(x)$ | error
---|---|---|---|---|---
10 | 4 | 3.21 | $-19.75\,\%$ | 3.08 | $-23.00\,\%$
$10^{2}$ | 26 | 22.67 | $-12.81\,\%$ | 25.58 | $-1.62\,\%$
$10^{3}$ | 201 | 185.11 | $-7.91\,\%$ | 202.61 | 1.80 %
$10^{4}$ | 1 680 | 1 603.11 | $-4.58\,\%$ | 1 710.35 | 1.80 %
$10^{5}$ | 14 902 | 14 338.63 | $-3.78\,\%$ | 15 069.0 | 1.12 %
$10^{6}$ | 135 124 | 130 893.21 | $-3.13\,\%$ | 136 274.75 | 0.85 %
$10^{7}$ | 1 243 370 | 1 211 835.68 | $-2.54\,\%$ | 1 253 639.87 | 0.83 %
$10^{8}$ | 11 587 149 | 11 335 684.78 | $-2.17\,\%$ | 11 672 710.45 | 0.74 %
$10^{9}$ | 108 941 388 | 106 873 861.02 | $-1.90\,\%$ | 109 666 579.94 | 0.67 %
$10^{10}$ | 1 031 330 156 | 1 013 894 469.43 | $-1.69\,\%$ | 1 037 530 754.16 | 0.60 %
Table 1. Estimates $cx/\sqrt{\ln x}$ and
$E(x)=c\int_{2}^{x}\frac{dt}{\sqrt{\ln t}}$ for $g(x)$.

An error of $1.69\%$ for $x=10^{10}$ is not impressive, but then the simple
form $\pi(x)\sim x/\ln x$ of the Prime Number Theorem, which can also be
viewed as an instance of Raikov’s Theorem, has an error of about $5\%$ here.
Motivated by the much greater accuracy of ${\rm Li}(x)=\int_{2}^{x}(\ln
t)^{-1}dt$ as an estimate for $\pi(x)$, we considered whether

$E(x):=c\negthickspace\int_{2}^{x}\negthickspace\frac{dt}{\sqrt{\ln t}}$

might give a better estimate for $g(x)$. The evidence is given in the fifth
and sixth columns of Table 1.

More generally, Raikov’s Theorem (7) can be restated as

$\sum_{n\leq x}a_{n}\sim
c\negthickspace\int_{2}^{x}\negthickspace\frac{dt}{(\ln t)^{\alpha}},$

where $c$ and $\alpha$ are as defined in Section 2, since the two estimates
are asymptotically equivalent. Our experience, together with heuristic
arguments based on the expected values of certain random variables, suggests
that this integral form will usually give more accurate approximations.

Digression. Motivated by curiosity rather than problems involving Riemann
surfaces, we note that if one defines

${\mathcal{P}}^{+}:=\\{q\mid q\;\hbox{is prime and}\;q\equiv 1\;{\rm
mod}\,(4)\\}$

then by essentially the same argument the corresponding counting function
$g^{+}(x)$ for $\mathcal{P}^{+}$-integers $n\leq x$ is estimated by a similar
formula

$g^{+}(x)\sim\frac{c^{+}x}{\sqrt{\ln
x}}\;=\sqrt{\frac{e^{\gamma}}{\pi}}\cdot\frac{C^{+}}{2}\cdot\frac{x}{\sqrt{\ln
x}}\quad\hbox{as}\quad x\to\infty,$

where

$C^{+}:=C(4,-1)=0.8689277682343238299091527791046529122939\ldots$

according to [10, 11]. Comparing this value with that of $C$ shows that
$\mathcal{P}$-integers appear nearly $50\%$ more frequently than
$\mathcal{P}^{+}$-integers. This is despite the fact that the numbers
$\pi^{\pm}(x)$ of primes $q\leq x$ satisfying $q\equiv\pm 1$ mod $(4)$ satisfy

$\pi^{+}(x)\sim\pi^{-}(x)\sim\frac{1}{2}\pi(x)\quad\hbox{as}\quad x\to\infty.$

The reason for this surprisingly large bias is that the smallest primes
$q\equiv 1$ mod $(4)$, such as $5,13,17,29,\ldots$, are almost always larger
than the corresponding primes congruent to $-1$ mod $(4)$, such as
$3,7,11,19,\ldots$ (see [5] for a very readable account of this and other
similar phenomena), so that $\mathcal{P}^{+}$-integers tend to be larger and
hence less frequently found than $\mathcal{P}$-integers. For example, the
first ten $\mathcal{P}$-integers are $3,7,9,11,19,21,23,27,31,33,$ while the
corresponding $\mathcal{P}^{+}$-integers are $5,13,17,25,29,37,41,53,61,65$.

Table 2 gives the numbers of $\mathcal{P}$\- and $\mathcal{P}^{+}$-integers
$n\leq x$ for various $x$. The constants $C(4,1)=1.29230415713\ldots$ and
$C(4,-1)=0.86892776823\ldots$ of interest to us have ratio
$C(4,-1)/C(4,1)=0.6723864220\ldots$. It looks plausible that the ratios in the
last column of the table tend to this value.

$x$ | $g(x)$ | $g^{+}(x)$ | $g^{+}(x)/g(x)$
---|---|---|---
10 | 4 | 2 | 0.5000000000
$10^{2}$ | 26 | 15 | 0.5769230769
$10^{3}$ | 201 | 123 | 0.6119402985
$10^{4}$ | 1 680 | 1 074 | 0.6392857142
$10^{5}$ | 14 902 | 9 623 | 0.6457522480
$10^{6}$ | 135 124 | 87 882 | 0.6503803913
$10^{7}$ | 1 243 370 | 814 183 | 0.6548195629
$10^{8}$ | 11 587 149 | 7 618 317 | 0.6574798511
$10^{9}$ | 108 941 388 | 71 838 469 | 0.6594231110
Table 2. Numbers $g(x)$ and $g^{+}(x)$ of $\mathcal{P}$\- and
$\mathcal{P}^{+}$-integers $n\leq x$.

It follows from a theorem of Iwaniec [7] that for each integer $A\neq 0$ the
number of primes $p\leq x$ of the form $\xi^{2}+\eta^{2}+A$ with
$\gcd(\xi,\eta)=1$ has order of magnitude $x/(\log x)^{3/2}$. Now a positive
integer $n$ is properly represented by the quadratic form $\xi^{2}+\eta^{2}$,
that is, with $\gcd(\xi,\eta)=1$, if and only if $n$ is not divisible by $4$
or by any prime $q\equiv-1$ mod $(4)$, or equivalently $n=m$ or $2m$ where $m$
is a ${\mathcal{P}}^{+}$-integer. Taking $A=-2$ we see that there are
infinitely many primes $p$ such that $p+2$ is a ${\mathcal{P}}^{+}$-integer.
Unfortunately, there is no quadratic form which plays a similar role for
$\mathcal{P}$-integers.

##  4\. An alternative evaluation of $c$

For an alternative approach to evaluating $c$ we will need the $L$-function

(16)
$L(s)=L(s,\chi):=\sum_{n=1}^{\infty}\frac{\chi(n)}{n^{s}}=\prod_{q}\left(1-\frac{\chi(q)}{q^{s}}\right)^{-1}$

corresponding to the non-principal character $\chi$ mod $(4)$ defined by
$\chi(n)=0$ for $n$ even, and $\chi(n)=\pm 1$ for $n\equiv\pm 1$ mod $(4)$.
According to [4, Exercise 11.8] this satisfies

(17) $\zeta(s)L(s)=\zeta_{K}(s),$

where $\zeta_{K}$ is the Dedekind zeta function of the field
$K={\mathbb{Q}}(i)$, defined by

(18)
$\zeta_{K}(s):=\sum_{I}\frac{1}{N(I)^{s}}=\prod_{P}\left(1-\frac{1}{N(P)^{s}}\right)^{-1}$

where $I$ and $P$ range over all ideals and prime ideals of
${\mathcal{O}}_{K}={\mathbb{Z}}[i]$, and $N(\;)$ denotes their norm. (To be
more specific, for $K={\mathbb{Q}}(i)$ the prime $q=2$ ramifies, primes
$q\equiv 1$ mod $(4)$ split into two prime ideals with norm $q$, and primes
$q\equiv-1$ mod $(4)$ are inert, each giving one prime ideal with norm
$q^{2}$, so

$\zeta_{K}(s)=\left(1-\frac{1}{2^{s}}\right)^{-1}\cdot\prod_{q\equiv
1(4)}\left(1-\frac{1}{q^{s}}\right)^{-2}\cdot\prod_{q\equiv-1(4)}\left(1-\frac{1}{q^{2s}}\right)^{-1};$

by using

$1-\frac{1}{q^{2s}}=\left(1-\frac{1}{q^{s}}\right)\left(1+\frac{1}{q^{s}}\right)$

one sees that this is equal to $\zeta(s)L(s)$.)

The Riemann zeta function $\zeta$ has a simple pole at $s=1$, with residue $1$
(see [4, Theorem 8.28], for example). By [4, Theorem 11.7] the function
$\zeta_{K}$, like $\zeta$, has a simple pole at $s=1$, but with residue
$\rho_{K}$ given by the class number formula

(19) $\rho_{K}=\lim_{s\to 1}(s-1)\zeta_{K}(s)=\frac{2\pi h}{w\sqrt{|D|}}$

where $D=-4$ is the discriminant of $K$, $h=1$ is its class number and $w=4$
is the number of units in ${\mathcal{O}}_{K}$, so that $\rho_{K}=\pi/4$.
(There is also a geometric proof of this, counting lattice points in circles
and squares: see [16] for an outline.) Since $\zeta$ and $\zeta_{K}$ have
meromorphic continuations in a neighbourhood of $1$, it follows that the
function $L=\zeta_{K}/\zeta$ is meromorphic there, with

(20) $L(1)=\lim_{s\to 1}\frac{\zeta_{K}(s)}{\zeta(s)}=\rho_{K}=\frac{\pi}{4}.$

We can now calculate $\alpha$. Comparing the Euler product expansions for
$\zeta$ and $L$ shows that

$\frac{\left(1-2^{-s}\right)\zeta(s)}{L(s)}=\prod_{q\equiv-1(4)}\frac{1+q^{-s}}{1-q^{-s}}=\prod_{q\equiv-1(4)}\frac{(1-q^{-2s})}{(1-q^{-s})^{2}}=\frac{F(s)^{2}}{F(2s)}.$

Since the functions $1-2^{-s}$, $L(s)$ and $F(2s)$ are analytic and non-zero
at $s=1$, and $\zeta(s)$ has a simple pole there, it follows that $F(s)^{2}$
also has a simple pole at $s=1$, so that $\alpha=\frac{1}{2}$. Thus

(21) $g(x)\sim\frac{cx}{\sqrt{\ln x}}\quad\hbox{as}\quad x\to\infty$

for some constant $c$.

We can now determine $c$. In Raikov’s Theorem we have

$c=\frac{H(1)}{\Gamma(1-\alpha)}=\frac{H(1)}{\Gamma(\frac{1}{2})}=\frac{H(1)}{\sqrt{\pi}},$

where $H(s)=\sqrt{s-1}F(s)$ near $s=1$, so that

$H(1)^{2}=\lim_{s\to 1}\,(s-1)F(s)^{2}=\frac{F(2)}{2L(1)}\lim_{s\to
1}\,(s-1)\zeta(s)=\frac{2F(2)}{\pi}$

since $L(1)=\pi/4$ and $\lim_{s\to 1}(s-1)\zeta(s)=1$. Thus

$c=\frac{H(1)}{\sqrt{\pi}}=\frac{1}{\pi}\sqrt{2F(2)}$

where

$F(2)=\prod_{q\equiv-1(4)}\left(1-\frac{1}{q^{2}}\right)^{-1}=\sum_{n=1}^{\infty}\frac{g_{n}}{n^{2}}=1+\frac{1}{3^{2}}+\frac{1}{7^{2}}+\cdots.$

The partial sums $\sum_{n=1}^{N}g_{n}/n^{2}$ of this series converge to $F(2)$
from below as $N\to\infty$. To obtain a sequence converging to $F(2)$ from
above we use

$\sum_{n=1}^{\infty}\frac{1}{n^{2}}=\zeta(2)=\frac{\pi^{2}}{6}$

(see [4, Exercise 9.7], for example), so

$\sum_{{\rm
odd}\,n}\frac{1}{n^{2}}=\zeta(2)-\frac{1}{2^{2}}\zeta(2)=\frac{\pi^{2}}{8}$

and hence

$F(2)=\frac{\pi^{2}}{8}-\sum_{n=1}^{\infty}\frac{h_{n}}{n^{2}}=\frac{\pi^{2}}{8}-\left(\frac{1}{5^{2}}+\frac{1}{13^{3}}+\cdots\right)$

where $h_{n}=1-g_{n}$ or $0$ as $n$ is odd or even, so that the sum is over
all odd $n\not\in{\mathcal{P}}^{*}$. Thus

(22)
$L:=\sum_{n=1}^{N}\frac{g_{n}}{n^{2}}<F(2)<\frac{\pi^{2}}{8}-\sum_{n=1}^{N}\frac{h_{n}}{n^{2}}\;:=U$

for all $N$, with both bounds converging monotonically to $F(2)$ as
$N\to\infty$, so one can evaluate $F(2)$ to any desired accuracy by taking $N$
sufficiently large.

For example, taking $N=10^{9}$ in (22) and evaluating terms to $30$ decimal
places gives upper and lower bounds

$U=1.16807558580668082351574608710,$ $L=1.16807558530668082351574605022,$

with $U-L=5.0000000000000003688\times 10^{-10}$. For the constant

$c=\frac{1}{\pi}\sqrt{2F(2)}$

this gives upper and lower bounds

$U^{\prime}=0.486519888468395460325765388132,$
$L^{\prime}=0.486519888364266948773906811661,$

with $U^{\prime}-L^{\prime}=1.04128511551858576471\times 10^{-10}$. The value

$c=0.4865198883858909971272456405868234055382\ldots$

given in (15) lies between these two bounds, and they agree with it in their
first nine significant figures.

Remark The values of $L$ and $U$ in (22) for any given $N$ depend on partial
sums of the Dirichlet series for the sequences $(g_{n})$ and $(h_{n})$,
representing odd integers $n$ which are or are not $\mathcal{P}$-integers. For
small $N$ the former predominate, so that the lower bounds $L$ are closer to
$F(2)$ than the upper bounds $U$ are; since both converge monotonically to
$F(2)$, this means that $U$ decreases faster than $L$ increases, especially as
non-$\mathcal{P}$-integers eventually start to predominate as $N$ increases.

To illustrate this, let us denote by $L_{k}$ and $U_{k}$ the lower and upper
bounds for $F(2)$ computed over the segment $n\leq N=10^{k}$, and let $R_{k}$
denote the ratio between the values of their rate of change, that is,

$R_{k}=\frac{L_{k}-L_{k-1}}{U_{k-1}-U_{k}}.$

For comparison, let $r_{k}$ denote the ratio $g(x)/h(x)$ of
$\mathcal{P}$-integers to non-$\mathcal{P}$-integers among the odd integers
$n\leq x=10^{k}$. These ratios $R_{k}$ and $r_{k}$ are shown, for
$k=1,\ldots,9$, in Table 3. It follows from Theorem 3.1 that $R_{k}$ and
$r_{k}$ converge to $0$ as $k\to\infty$, and Table 3 shows that they do so ‘in
step’.

$k$ | $R_{k}$ | $r_{k}$
---|---|---
1 | 3.5966238347 | 4.0000000000
2 | 0.9856093680 | 1.0833333333
3 | 0.6579817478 | 0.6722408026
4 | 0.5213884648 | 0.5060240963
5 | 0.4445379068 | 0.4245825972
6 | 0.3809371214 | 0.3703285499
7 | 0.3382979524 | 0.3309801604
8 | 0.3076691054 | 0.3016477220
9 | 0.2834484269 | 0.2785807156
Table 3. Comparison of ratios $R_{k}$ and $r_{k}$.

##  5\. The twin primes conjecture

It is useful now to recall the heuristic arguments used by Hardy and
Littlewood to justify their work in [6] on, among various other problems, the
Twin Primes Conjecture, which asserts that there are infinitely many pairs of
twin primes $p,p+2$.

By the Prime Number Theorem, the probability of a natural number $p\leq x$
being prime is asymptotic to $x/\ln x$, so the numbers of primes and of pairs
of primes $p,p^{\prime}\leq x$ are asymptotically equivalent to

${\rm Li}(x):=\int_{2}^{x}\negthickspace\frac{dt}{\ln t}\quad\hbox{and}\quad
I_{2}(x):=\int_{2}^{x}\negthickspace\frac{dt}{(\ln t)^{2}}.$

However, if $p$ and $p^{\prime}$ are related by an equation such as
$p^{\prime}=p+2$, rather than chosen independently, the estimate $I_{2}(x)$
will be incorrect. Since natural numbers are uniquely determined by their
residues mod $(q)$ for all primes $q$ it makes sense to apply the restriction
$p^{\prime}=p+2$ through its effect on these residues. For each prime $q$, we
need to consider the probabilities that $p$ and $p^{\prime}$, chosen
independently or with $p^{\prime}=p+2$, being both prime, are both coprime to
$q$ (we neglect the vanishingly small probability that either of them is equal
to $q$). By a simple count of congruence classes mod $(q)$ we see that these
probabilities are respectively

$\left(1-\frac{1}{q}\right)^{2}\quad\hbox{and}\quad
1-\frac{\omega_{f}(q)}{q},$

where $\omega_{f}(q)$ is the number of roots of the polynomial $f(t)=t(t+2)$
mod $(q)$. In order to replace the first probability, implicit in the estimate
$I_{2}(x)$, with the second more correct probability, we therefore multiply
$I_{2}(x)$ by the correction factor

$C_{2}(q)=\left(1-\frac{1}{q}\right)^{-2}\left(1-\frac{\omega_{f}(q)}{q}\right).$

Doing this for each prime $q$ we multiply $I_{2}(x)$ by a correction term,
called a Hardy–Littlewood constant,

$C(f)=\prod_{q\;{\rm prime}}C_{2}(q)=\prod_{q\;{\rm
prime}}\left(1-\frac{1}{q}\right)^{-2}\left(1-\frac{\omega_{f}(q)}{q}\right).$

Clearly $\omega_{f}(2)=1$ while $\omega_{f}(q)=2$ for each prime $q>2$, so we
have an estimate for the number $\pi_{2}(x)$ of twin prime pairs $p,p+2\leq x$
of the form

(23)
$E_{2}(x)=C(f)I_{2}(x)=2C_{2}I_{2}(x)=2C_{2}\negthickspace\int_{2}^{x}\negthickspace\frac{dt}{(\ln
t)^{2}}\,,$

where the initial factor $2$ is the correction factor $C_{2}(2)$ corresponding
to the prime $q=2$, and

(24)
$C_{2}:=\prod_{q>2}\left(1-\frac{1}{q}\right)^{-2}\left(1-\frac{2}{q}\right)=\prod_{q>2}\left(1-\frac{1}{(q-1)^{2}}\right)=0.6601618158\ldots$

with the product over all odd primes $q$.

The conjecture is that

$\pi_{2}(x)\sim E_{2}(x)\quad\hbox{as}\quad x\to\infty.$

Although it is unproved, there is strong evidence for this conjecture. For
example, it is known that

$\pi_{2}(10^{18})=808\,675\,888\,577\,436,$

while Maple evaluates

$E_{2}(10^{18})=808\,675\,901\,493\,606.3\ldots.$

The relative error is $0.0000016\%$.

In 1962 Bateman and Horn [1] made a wide-ranging generalisation of the
Hardy–Littlewood estimate, giving a similar conjectured estimate $E(x)$ for
the number $Q(x)$ of natural numbers $t\leq x$ at which a given finite set of
polynomials $f_{1}(t),\ldots,f_{k}(t)\in{\mathbb{Z}}[t]$ simultaneously take
prime values. It assumes that these polynomials satisfy the following
conditions, which are obviously necessary for there to be infinitely many such
$t$ (Schinzel’s Hypothesis [15] asserts that they are also sufficient):

  1. (1)

each $f_{i}(t)$ has a positive leading coefficient;

  2. (2)

each $f_{i}(t)$ is irreducible in ${\mathbb{Z}}[t]$;

  3. (3)

the product $f(t):=f_{1}(t)\ldots f_{k}(t)$ is not identically zero modulo any
prime.

The Bateman–Horn Conjecture (BHC) asserts that

(25) $Q(x)\sim E(x):=C\negthickspace\int_{a}^{x}\negthickspace\frac{dt}{\ln
f_{1}(t)\ldots\ln f_{k}(t)}\quad\hbox{as}\quad x\to\infty,$

where $C$ is the Hardy–Littlewood constant

(26) $C=C(f_{1},\ldots,f_{k})=\prod_{q\;{\rm
prime}}\left(1-\frac{1}{q}\right)^{-k}\left(1-\frac{\omega_{f}(q)}{q}\right),$

$\omega_{f}(q)$ is the number of roots of $f$ mod $(q)$, and the lower limit
$a$ in the integral in (25) is chosen so that the integral avoids
singularities where $\ln f_{i}(t)=0$ for some $i$. If conditions (1) to (3)
are satisfied then the infinite product in (26) converges to a limit $C>0$,
while the definite integral in (25) diverges to $+\infty$ with $x$, so if the
BHC is true then $Q(x)\to+\infty$ and Schinzel’s Hypothesis, that the
polynomials $f_{i}$ simultaneously take prime values for infinitely many
$t\in\mathbb{N}$, is verified.

This includes the cases of the twin primes and the Sophie German primes
conjectures, where $f_{1}(t)=t$ and $f_{2}(t)=t+2$ or $2t+1$ respectively. The
BHC has been proved only in the case of a single polynomial of degree 1: this
is the quantified version, due to de la Vallée Poussin, of Dirichlet’s Theorem
on primes in an arithmetic progression $at+b$. Nevertheless, the estimates
produced by the BHC, in a wide range of applications, agree remarkably well
with observed counts obtained by primality-testing. As an example, with
$f_{1}(t)=4t+1$ and $f_{2}(t)=4t+3$, Table 4 shows the BHC estimates
$E_{2}^{+}(x)$ for the numbers $\pi_{2}^{+}(x)$ of twin primes $p,p+2\leq x$
with $p\equiv 1$ mod $(4)$, those yielding BCC primes $p$.

$x$ | $\pi_{2}^{+}(x)$ | $E_{2}^{+}(x)$ | error
---|---|---|---
10 | 1 | 1.148985018 | 14.8985 %
$10^{2}$ | 4 | 5.498634634 | 37.4659 %
$10^{3}$ | 19 | 21.62864106 | 13.8350 %
$10^{4}$ | 105 | 105.8363607 | 0.79653 %
$10^{5}$ | 604 | 623.0852586 | 3.15981 %
$10^{6}$ | 4 046 | 4 122.745734 | 1.89683 %
$10^{7}$ | 29 482 | 29 375.63915 | $-0.3608\,\%$
$10^{8}$ | 220 419 | 220 182.6280 | $-0.1072\,\%$
$10^{9}$ | 1 712 731 | 1 712 652.809 | $-0.0046\,\%$
$10^{10}$ | 13 706 592 | 13 705 706.99 | $-0.0065\,\%$

Table 4. Numbers of pairs of twin primes $p,p+2\leq x$ with $p\equiv 1$ mod
$(4)$.

”

##  6\. In the footsteps of Hardy and Littlewood

We will now adapt the heuristic arguments used by Hardy and Littlewood, and
later by Bateman and Horn, to the slightly different context of our BCC primes
problem.

Using the Prime Number Theorem and Theorem 3.1 to give the distributions of
prime numbers $p$ and of $\mathcal{P}$-integers $n$, we first consider

$I(x):=\int_{2}^{x}\negthickspace\frac{c\,dt}{(\ln t)^{3/2}}$

as an estimate for the number $a(x)$ of BCC primes $p\leq x$, where $c$ is as
in (15). Of course, this treats $p$ and $n$ as independent random variables,
and takes no account of the fact that $n=p+2$. In following the example of
Hardy and Littlewood, with a similar heuristic justification, we now apply
correction factors only for primes $q\not\in{\mathcal{P}}$, that is, $q=2$ or
$q\equiv 1$ mod $(4)$, since it is for such $q$ that we need to replace the
probability that independent variables $p$ and $n$ are both coprime to $q$
with the corresponding probability for $p$ and $p+2$. We therefore multiply
this estimate $I(x)$ by $2C_{2}^{+}$ where

(27) $C_{2}^{+}:=\prod_{q\equiv
1(4)}\left(1-\frac{1}{q}\right)^{-2}\left(1-\frac{2}{q}\right)=\prod_{q\equiv
1(4)}\left(1-\frac{1}{(q-1)^{2}}\right),$

with the extra factor $2$ corresponding to the prime $q=2$.

This leads to the following conjecture:

######  Conjecture 6.1.

The function $a(x)$ which counts BCC primes $p\leq x$ satisfies

(28) $a(x)\sim
2C_{2}^{+}I(x)=c^{\prime}\negthickspace\int_{2}^{x}\negthickspace\frac{dt}{(\ln
t)^{3/2}}\quad\hbox{as}\quad x\to\infty,$

where

(29) $c^{\prime}=2C_{2}^{+}c=\prod_{q\equiv
1\,(4)}\left(1-\frac{1}{(q-1)^{2}}\right)\cdot\sqrt{\frac{e^{\gamma}}{\pi}}\cdot
C(4,1).$

##  7\. Computations

The various estimates we have discussed can be tested computationally by using
Maple. The definite integrals, such as those appearing in (23) and (28), are
evaluated accurately and rapidly, even for large values of $x$, by numerical
integration. The Hardy–Littlewood constants, such as those in (24) and (27),
are defined as infinite products which converge conditionally and slowly; good
approximations can be found by taking partial products over large initial
segments of the relevant primes, in increasing order, but on a laptop these
calculations can take a matter of hours. Having evaluated the various
estimates, one can compare them with the actual numbers of terms being counted
by using the Rabin–Miller primality test within Maple.

Using Maple to evaluate (27) for primes $q\leq 10^{8}$ we find that

$C_{2}^{+}=0.9230611326\ldots$

The values of $\sqrt{e^{\gamma}/\pi}$ and $C(4,1)$ given in Section 3 allow us
to deduce that

(30) $c^{\prime}=0.8981751984\ldots$

To test the value for $c^{\prime}$ in (30) let us define $C(x)$ by

$a(x)=C(x)\negthickspace\int_{2}^{x}\negthickspace(\ln t)^{-3/2}dt,$

so we expect that $C(x)\to c^{\prime}$ as $x\to\infty$. Using Maple to
evaluate $a(x)$ and $\int_{2}^{x}(\ln t)^{-3/2}dt$, we found the values of
$C(x)$ shown in Tables 5 and 6. We note that, after some initial instability
due to the relatively small numbers of BCC primes appearing, the values of
$C(x)$ decrease almost monotonically towards $c^{\prime}$ from $x=10^{7}$ to
$x=35\cdot 10^{9}$.

$x$ | $a(x)$ | $C(x)$
---|---|---
$10^{1}$ | 2 | 0.4688555840
$10^{2}$ | 12 | 0.7153107401
$10^{3}$ | 65 | 0.8472363117
$10^{4}$ | 388 | 0.8706477077
$10^{5}$ | 2 708 | 0.9004062930
$10^{6}$ | 19 969 | 0.9024565742
$10^{7}$ | 155 369 | 0.9040719795
$10^{8}$ | 1 250 182 | 0.9023589943
$10^{9}$ | 10 345 920 | 0.9011815839
$10^{10}$ | 87 545 946 | 0.9010054301
Table 5. Evolution of the coefficients $C(x)$ for $x=10^{k}$, $k=1,\ldots,10$. $x$ | $a(x)$ | $C(x)$
---|---|---
$11\cdot 10^{9}$ | 95 675 252 | 0.9010083668
$12\cdot 10^{9}$ | 103 758 501 | 0.9010266288
$13\cdot 10^{9}$ | 111 794 166 | 0.9010090565
$14\cdot 10^{9}$ | 119 795 477 | 0.9010345605
$15\cdot 10^{9}$ | 127 757 388 | 0.9010373485
$16\cdot 10^{9}$ | 135 683 004 | 0.9010241227
$17\cdot 10^{9}$ | 143 578 133 | 0.9010192899
$18\cdot 10^{9}$ | 151 444 525 | 0.9010197264
$19\cdot 10^{9}$ | 159 283 669 | 0.9010228308
$20\cdot 10^{9}$ | 167 092 278 | 0.9010018195
$21\cdot 10^{9}$ | 174 878 590 | 0.9009949342
$22\cdot 10^{9}$ | 182 641 563 | 0.9009888293
$23\cdot 10^{9}$ | 190 382 424 | 0.9009836658
$24\cdot 10^{9}$ | 198 100 186 | 0.9009699991
$25\cdot 10^{9}$ | 205 796 174 | 0.9009503295
$26\cdot 10^{9}$ | 213 474 480 | 0.9009388161
$27\cdot 10^{9}$ | 221 136 287 | 0.9009360946
$28\cdot 10^{9}$ | 228 778 913 | 0.9009277506
$29\cdot 10^{9}$ | 236 403 538 | 0.9009161678
$30\cdot 10^{9}$ | 244 014 302 | 0.9009145456
$31\cdot 10^{9}$ | 251 607 900 | 0.9009079931
$32\cdot 10^{9}$ | 259 184 139 | 0.9008943241
$33\cdot 10^{9}$ | 266 743 992 | 0.9008756612
$34\cdot 10^{9}$ | 274 293 778 | 0.9008715586
$35\cdot 10^{9}$ | 281 827 468 | 0.9008601041
Table 6. Evolution of the coefficients $C(x)$ for $x=m\cdot 10^{9}$,
$m=11,\ldots,35$.

Beyond this point, in view of the modest computing facilities available we
considered segments $[x,y]\subset\mathbb{R}$ of length $y-x=10^{8}$ or
$10^{7}$, starting at values $x=10^{15},10^{20},\ldots,10^{50}$. The results
are shown in Table 7, where $a(x,y)=a(y)-a(x)$ is the number of BCC primes
$p\in[x,y]$, and $C(x,y)$ is defined by the equation

$a(x,y)=C(x,y)\cdot\negthinspace\int_{x}^{y}\negthickspace\frac{dt}{(\ln
t)^{3/2}}.$

Initially we see a further monotonic decrease towards $c^{\prime}$, but then
$C(x,y)<c^{\prime}$ when $x=10^{30}$ and $y=10^{30}+10^{8}$. We are not very
concerned about this: for instance, Littlewood [12] famously showed that the
error ${\rm Li}(x)-\pi(x)$ in the Prime Number Theorem changes sign infinitely
many times (see also [14] for a similar phenomenon related to Mertens’s Third
Theorem), so why not here? The instability in the last four rows is perhaps
due to the relatively small numbers of BCC primes appearing in these shorter
intervals: for instance, compare the values of $a(x,y)$ with that of $a(x)$
for the equivalent interval $[0,10^{7}]$.

$x$ | $y-x$ | $a(x,y)$ | $C(x,y)$
---|---|---|---
$10^{15}$ | $10^{8}$ | 442 649 | 0.8985037831
$10^{20}$ | $10^{8}$ | 287 429 | 0.8982538979
$10^{25}$ | $10^{8}$ | 205 949 | 0.8994835673
$10^{30}$ | $10^{8}$ | 156 398 | 0.8979178604
$10^{35}$ | $10^{7}$ | 12 389 | 0.8963174519
$10^{40}$ | $10^{7}$ | 10 299 | 0.9103503807
$10^{45}$ | $10^{7}$ | 8 554 | 0.9022181293
$10^{50}$ | $10^{7}$ | 7 507 | 0.9273527313
Table 7. Estimation over segments $[x,y]$.

##  8\. Some remarks

1. The estimate $g(x)\sim cx/(\ln x)^{1/2}$ in Theorem 3.1, including the evaluation of $c$ in (15), is proved, as a consequence of Raikov’s Theorem. The estimate $a(x)\sim c^{\prime}\int_{2}^{x}(\ln t)^{-3/2}dt$ in Conjecture 6.1, including the formula for $c^{\prime}$ given there, is just a conjecture, and we do not expect to see a proof for it soon. Nevertheless, we feel that the computational data in Section 7 give plausible evidence for the validity of this estimate, and hence for the infinitude of BCC primes.

2. In the construction of Riemann surfaces used in [2] and considered here, it is important for group-theoretic reasons that the prime $p$ and the $\mathcal{P}$-integer $n$ should differ by $2$. However, it is clear that in the arguments we have used one could generalise this relationship, and still obtain similar estimates for the distribution of such primes. For instance, one could replace the difference $2$ here with any non-zero even integer. The following section gives a less trivial example, also arising from [2].

##  9\. A similar problem

In [2], having obtained an Accola–Maclachlan-type lower bound $M(g)\geq
4(g-1)$ for the largest possible order $M(g)$ of the automorphism group of a
pseudo-real Riemann surface of odd genus $g\geq 3$, Bujalance, Cirre and
Conder have presented a set $\mathcal{G}$ of genera $g$ for which this bound
is sharp. (The sharpness of the corresponding bound $M(g)\geq 2g$ for even
genera $g\geq 2$ is left as an open problem.) For this they give a similar
construction (see the remark following their Theorem 5.2) of a family of
pseudo-real Riemann surfaces of genus $g=2p+1$ where $p$ is what we will call
a BCC2 prime, defined as follows:

######  Definition 9.1.

A BCC2 prime is a prime number $p$ such that

(31) $p\equiv 3\;{\rm mod}\,(8),\;p\equiv 2\;\hbox{or}\;5\;{\rm
mod}\,(9)\;\hbox{and}\;p\not\equiv 5\;{\rm mod}\,(7),$

and

(32) $n:=p+1\;\hbox{is not divisible by}\;11,23,47\;\hbox{or any
prime}\;q\equiv 1\;\hbox{mod}\;(3)\;\hbox{or}\;q\equiv 1\;\hbox{mod}\;(4).$

This raises the following problem, which appears to be even more challenging
than Problem 1.2:

######  Problem 9.2.

Are there infinitely many BCC2 primes?

These primes are much rarer than the BCC primes considered earlier: for
example, compare Table 8, which shows the first sixteen of them, with the
corresponding list of BCC primes in (1).

$p$ | factors of $n=p+1=12m$ |  $m$ mod $(7)$ |  $r=m$ mod $(84)$
---|---|---|---
$11$ | $2^{2}\cdot 3$ |  |
$1283$ | $2^{2}\cdot 3\cdot 107$ | $2$ | $23$
$1571$ | $2^{2}\cdot 3\cdot 131$ | $5$ | $47$
$2003$ | $2^{2}\cdot 3\cdot 167$ | $6$ | $83$
$3011$ | $2^{2}\cdot 3\cdot 251$ | $6$ | $83$
$7043$ | $2^{2}\cdot 3\cdot 587$ | $6$ | $83$
$7907$ | $2^{2}\cdot 3\cdot 659$ | $1$ | $71$
$8627$ | $2^{2}\cdot 3\cdot 719$ | $5$ | $47$
$9923$ | $2^{2}\cdot 3\cdot 827$ | $1$ | $71$
$10\,067$ | $2^{2}\cdot 3\cdot 839$ | $6$ | $83$
$15\,107$ | $2^{2}\cdot 3\cdot 1259$ | $6$ | $83$
$15\,683$ | $2^{2}\cdot 3\cdot 1307$ | $6$ | $83$
$17\,123$ | $2^{2}\cdot 3\cdot 1427$ | $6$ | $83$
$17\,987$ | $2^{2}\cdot 3\cdot 1499$ | $1$ | $71$
$18\,131$ | $2^{2}\cdot 3\cdot 1511$ | $6$ | $83$
$19\,427$ | $2^{2}\cdot 3\cdot 1619$ | $2$ | $23$
Table 8. The first sixteen BCC2 primes $p$.

First let us restate the definition of BCC2 primes in terms of $n$. The
congruences mod $(8)$ and mod $(9)$ in condition (31) are equivalent to
$n=p+1\equiv 4$ mod $(8)$ and $n\equiv 3$ or $6$ mod $(9)$, that is, $n=12m$
for some $m$ coprime to $12$; when this is satisfied the condition
$p\not\equiv 5$ mod $(7)$ in (31) is equivalent to $m\not\equiv 4$ mod $(7)$.
Thus condition (31) is equivalent to

(33) $n=12m\;\hbox{where}\;(12,m)=1\;\hbox{and}\;m\not\equiv 4\,{\rm
mod}\,(7).$

Similarly, when this is satisfied, condition (32) is equivalent to

(34) $\hbox{each prime $q$ dividing $m$ satisfies}\;47<q\equiv-1\,{\rm
mod}~{}(12).$

Thus conditions (31) and (32), taken together, are equivalent to (33) and
(34), also taken together. Recall that in addition we require $n-1$ ($=p$) to
be prime (and greater than $7$), which implies that $m\not\equiv 3$ mod $(7)$
since otherwise we would have $12m-1\equiv 36-1\equiv 0$ mod $(7)$, while
conditions (33) and (34) imply that $m\not\equiv 4$ or $0$ mod $(7)$, so we
have

(35) $m\equiv 1,2,5\;\hbox{or}\;6\;\hbox{mod}\;(7).$

In Table 8, apart from $p=11$, all entries have the property that $n=12m$ for
some prime $m$. Indeed, the first exception to this is $p=41\,771$, with
$n=12\cdot 59^{2}$. Moreover, among the $172\,515$ BCC2 primes $p\leq 10^{9}$,
there are only $41\,711$ for which $m$ is not prime. Motivated by this,
instead of attempting to estimate the distribution of the whole set of BCC2
primes, as we did in the case of the BCC primes, we will apply the
Bateman–Horn Conjecture (see Section 5) to give strong evidence that there are
infinitely many of them for which $m$ is prime. (This is analogous to the
argument at the end of Section 5, where we considered those BCC primes $p$
which are members of twin primes $p,p+2$.)

Let us therefore assume that $n=12m$ for some prime $m>47$ such that
$m\equiv-1$ mod $(12)$ and $m\equiv 1,2,5$ or $6$ mod $(7)$, or equivalently
$m\equiv r$ mod $(84)$ where $r=71,23,47$ or $83$ by the Chinese Remainder
Theorem. Clearly, any such choice of $m$ satisfies conditions (33), (34) and
(35). For each such $r$ we are therefore looking for integers $t\geq 0$ such
that the polynomials

$f_{1}(t)=84t+r$

and

$f_{2}(t)=12f_{1}(t)-1=1008t+12r-1$

both take prime values, namely $m$ and $p$. These two polynomials are
irreducible, with positive leading coefficients, and their product $f$ is not
identically zero modulo any prime, so we have four instances of the
Bateman–Horn Conjecture, one for each value of $r$. In each case the BHC gives
an asymptotic estimate $E_{r}(x)$ for the number $Q_{r}(x)$ of integers $t\leq
x$ such that $f_{1}(t)$ and $f_{2}(t)$ are both prime. We have
$\omega_{f}(q)=0$ for the primes $q=2,3$ and $7$, and $\omega_{f}(q)=2$ for
all other primes, so the Hardy–Littlewood constant $C=C(f_{1},f_{2})$, which
is independent of $r$, is positive. In fact, comparing the infinite products
for $C$ and for $C_{\rm twins}=2C_{2}$ (see Section 5), which differ at only
these three primes, shows that

$C=\frac{42}{5}\cdot C_{\rm twins}=\frac{42}{5}\cdot
1.3203236316\ldots=11.0907185062\ldots.$

It follows that $E_{r}(x)\to+\infty$ as $x\to\infty$, giving evidence that for
each $r$ there are infinitely many prime values $m$ and $p$ for $f_{1}$ and
$f_{2}$, and thus infinitely many BCC2 primes $p$.

To test this approach we found that of the $172\,515$ BCC2 primes $p\leq
10^{9}$, the number with $m$ prime, as we have been assuming, is $Q=130\,804$.
The four BHC estimates

$E_{r}(x)=C\cdot\negthinspace\int_{2}^{x}\negthickspace\frac{dt}{\ln
f_{1}(t)\ln f_{2}(t)}$

were computed by taking $x=992\,062$; beyond this point the primes
$p=1008t+12r-1$ begin to be greater than $10^{9}$. The four values obtained
were almost identical, ranging from $32\,709.13$ to $32\,709.29$. Their sum
gave an estimate of $E=130\,836.81$ for $Q$, with an error of $0.025\%$.

##  10\. Appendix: the pseudo-real surfaces

Here we briefly summarise the construction of the pseudo-real Riemann surfaces
in [2] which give rise to Problems 1.2 and 9.2. In the first case, for any
even $g\geq 2$ let $\Gamma$ be an NEC group with signature
$(1;-;[2,2,g];\\{-\\})$. This means that $\Gamma$ acts as a group of
isometries of the hyperbolic plane $\mathbb{H}$, with canonical generators $d$
(a glide reflection) and $x_{1},x_{2},x_{3}$ (elliptic elements), and defining
relations

$x_{1}^{2}=x_{2}^{2}=x_{3}^{g}=d^{2}x_{1}x_{2}x_{3}=1,$

so that the quotient-surface ${\mathbb{H}}/\Gamma$ is the real projective
plane with three cone-points of orders $2,2$ and $g$.

Now let $\theta:\Gamma\to{\rm C}_{2g}=\langle u\mid u^{2g}=1\rangle$ be the
surface-kernel epimorphism defined by

$d\mapsto u,\;x_{1}\mapsto u^{g},\;x_{2}\mapsto u^{g},\;x_{3}\mapsto u^{-2},$

and let $K=\ker\,\theta$. Then $K$ is contained in the orientation-preserving
subgroup $\Gamma^{+}$ of index $2$ in $\Gamma$, so $S:={\mathbb{H}}/K$ is a
compact Riemann surface, of genus $g$ by the Riemann–Hurwitz formula. There is
a natural action of $\Gamma/K$ and hence of ${\rm C}_{2g}$ on $S$, with the
elements of the subgroup $\langle u^{2}\rangle\cong{\rm C}_{g}$ and its
complement acting as conformal and anticonformal automorphisms of $S$, and
with $S/{\rm C}_{2g}\cong{\mathbb{H}}/\Gamma$; the canonical double cover
$S/{\rm C}_{g}\cong{\mathbb{H}}/\Gamma^{+}$ of this orbifold is the Riemann
sphere with four cone-points of order $2$ and two of order $g$.

The generator $u$ of ${\rm C}_{2g}$ acts as an anticonformal automorphism of
$S$ of order $2g$, whereas the only involution in $G$ (namely $u^{g}$) acts
conformally. One can choose the cone-points so that $\Gamma$ is a maximal NEC
group, in which case $S$ has no further automorphisms (conformal or
anticonformal), and is therefore pseudo-real. A detailed argument in [2,
Theorem 4.8] shows that provided $g$ satisfies various other conditions
(namely those defining the set $\mathcal{A}$, together with $g>30$), no
pseudo-real surface of genus $g$ can have an abelian group of conformal
automorphisms with more than $g$ elements, so that $M^{+}_{\rm ab}(g)=g$.

The construction of the surfaces $S$ giving rise to Problem 9.2 is similar,
except that in this case $\Gamma$ has signature $(1;-;[2,2,2];\\{-\\})$, $S$
has genus $g=2p+1$ with a full automorphism group

$\langle u,v\mid u^{2}=v^{4p}=1,v^{u}=v^{2p-1}\rangle\cong{\rm
C}_{4p}\rtimes{\rm C}_{2}$

of order $8p=4(g-1)$, and $\theta$ is given by

$d\mapsto uv,\;x_{1}\mapsto u,\;x_{2}\mapsto u,\;x_{3}\mapsto v^{2p}.$

It is shown in [2, Theorem 5.3] that if $p$ is a prime satisfying the
conditions stated in Section 9 then no pseudo-real surface of genus $g$ has a
larger automorphism group, so that $M(g)=4(g-1)$.

Acknowledgments The authors are grateful to Emilio Bujalance, Javier Cirre and
Marston Conder for creating the problems addressed in this paper, and to Ken
Ford for some very helpful advice.

## References

  * [1] P. T. Bateman and R. A. Horn, A heuristic asymptotic formula concerning the distribution of prime numbers, Math. Comp. 16 (1962), 220–228.
  * [2] E. Bujalance, F. J. Cirre and M. D. E. Conder, Bounds on the orders of groups of automorphisms of a pseudo-real surface of given genus, J. London Math. Soc. (2) 101 (2020), 877–906.
  * [3] A. C. Cojocaru and M. R. Murty, Introduction to Sieve Methods and their Applications, London Math. Soc. Student Texts, Cambridge University Press, 2005.
  * [4] G. Everest and T. Ward, Introduction to Number Theory, Graduate Texts in Math. 232, Springer, 2005.
  * [5] A. Granville and G. Martin, Prime number races, Amer. Math. Mothly 113 (2006), 1–33; see also arXiv: 0408319.
  * [6] G. H. Hardy and J. E. Littlewood, Some problems of ‘Partitio numerorum’; III: On the expression of a number as a sum of primes, Acta Math. 114 (1923), 215–273.
  * [7] H. Iwaniec, Primes of the type $\varphi(x,y)+A$, where $\varphi$ is a quadratic form, Acta Arithmetica 21 (1972), 203–234.
  * [8] E. Landau, Lösung des Lehmer’schen Problems, Amer. J. Math. 31 (1909), 86–102.
  * [9] A. Languasco and A. Zaccagnini, On the constant in the Mertens product for arithmetic progressions, I. Identities, Functi. Approx. Comment. Math. 42 (2010), 17–27; see also arXiv:0706.2807v4.
  * [10] A. Languasco and A. Zaccagnini, On the constant in the Mertens product for arithmetic progressions, II. Numerical values, Math. Comp. 78 (2009), 315–326; see also arXiv:0712.0665v1.
  * [11] A. Languasco and A. Zaccagnini, Computation of the Mertens constants in arithmetic progressions, www.math.unipd.it/ languasc/MCcomput.html.
  * [12] J. E. Littlewood, Sur la distribution des nombres premiers, Comptes Rendus 158 (1914), 1869–1872.
  * [13] D. A. Raikov, Generalisation of a theorem of Ikehara–Landau (in Russian), Mat. Sbornik 45 (1938), 559–68.
  * [14] G. Robin, Sur l’ordre maximum de la fonction somme des diviseurs, Séminaire Delange–Pisot–Poitou, Théorie des Nombres (1981–1982), Progress in Mathematics 38 (1983), 233–244.
  * [15] A. Schinzel and W. Sierpiński, Sur certaines hypothèses concernant les nombres premiers, Acta Arithmetica 4 (1958), no. 2, 185–208.
  * [16] Wikipedia “Class number formula”.
  * [17] Wikipedia “Gauss’s constant”.


# Convolutional neural network for self-mixing interferometric displacement
sensing

Stéphane Barland 1,* and François Gustave2  1Université Côte d’Azur - CNRS,
Institut de Physique de Nice, 1361 route des Lucioles, F-06560, Valbonne,
France
2DOTA, ONERA, Université Paris-Saclay, F-91123, Palaiseau, France
<EMAIL_ADDRESS>
††journal: oe

Self mixing interferometry is a well established interferometric measurement
technique. In spite of the robustness and simplicity of the concept,
interpreting the self-mixing signal is often complicated in practice, which is
detrimental to measurement availability. Here we discuss the use of a
convolutional neural network to reconstruct the displacement of a target from
the self mixing signal in a semiconductor laser. The network, once trained on
periodic displacement patterns, can reconstruct arbitrarily complex
displacement in different alignment conditions and setups. The approach
validated here is amenable to generalization to modulated schemes or even to
totally different self mixing sensing tasks.

##  1 Introduction

Optical interferometric measurements are routinely used in science and
engineering and many schemes can be used to adapt the approach to the specific
measurement to be performed. One particularly interesting and well established
method is the so-called self-mixing interferometry, which consists in
realizing interference between the beam reflected by a target and a reference
beam inside the laser resonator emitting the reference beam (see eg [1, 2, 3,
4, 5] for reviews). For its simplicity and versatility, many applications have
been envisioned and perhaps the most immediate is that of displacement
measurement. Two limit regimes are considered [6]: that of very small
displacement (much smaller than the laser wavelength) or the opposite case
where the displacement takes place over a very large number of wavelengths. In
the first case, information about the target displacement can be retrieved
from fitting the shape of the interferometric signal. In the latter case, most
of the information is obtained by counting the fringes that are observed as a
sawtooth signal whose symmetry depends on the direction of the motion. Despite
its apparent simplicity, this analysis is often complicated since the exact
shape of the interferometric signal depends on many factors including bias
current, target reflectivity, alignment conditions [7], and modal structure of
the laser which may even lead to a double-peak structure in each fringe [8].
Furthermore, on diffusive targets, speckle leads to an effective variation of
the feebdack parameters and therefore a change of the signal shape in the
course of the measurement. In practice, all these effects tremendously affect
the availability of self-mixing measurement setups. This has led to a number
of hardware and software proposals to either improve the signal quality or the
retrieval of the displacement from the interferometric signal [9, 10, 11, 12,
13, 14].

Computer neural networks are one of the many architectures which can be used
for machine learning tasks, whereby a computer is used to infer rules from a
set of data and results instead of providing results on the basis of an input
and a priori known rules. The training of a neural network leads to the
formulation of kind of statistical model [15], able to predict new results on
the basis of new data. These neural networks are already very widely used in
everyday life and they are proving increasingly useful many areas of research
and technology. In the specific context of interferometry, very few attempts
exist to date. They have been used to identify and count fringes in [16, 17,
18]. In [19] and [20] they have been used to pre-process self mixing traces
and in [21] they are used as a part of a self-mixing blood pressure
measurement scheme.

In the following, we discuss the use of a convolutional neural network for the
direct reconstruction of a displacement signal across many different alignment
conditions. We address a particularly delicate regime which is the one of "few
wavelengths" displacement. The neural network is first trained on a set of
periodic data and in different alignment conditions. Its performance in
reconstructing the displacement of a target from a self-mixing interferometric
signal is then validated on aperiodic times series whose continuous spectrum
spans more than three octaves and under different alignment conditions, not
used during training. We also analyze the robustness of the reconstruction to
the presence of very strong detection noise. Finally, we observe that the
neural network can, without any tuning, provide a sensible reconstruction of
the displacement of a target obtained on a different experimental setup based
on the same operating principle. We then briefly discuss some details on the
operation of such a network and some further possible uses of this approach in
the context of self mixing. Thus, a reasonably simple neural network such as
the one used here can become one of the tools which contribute to the
robustness and high availability of self-mixing interferometric setups.

##  2 Neural network design and training

###  2.1 Experimental setup

Figure 1: The self-mixing experimental setup consists of a laser diode (LD),
short focal length collimator (Col.), laser diode current driver (LD driver)
and voltage amplifier (RF amp).

The experimental arrangement, presented in Fig. 1, consists of a single
transverse mode laser emitting at $\lambda=1310~{}nm$ (ML725B8F) whose
threshold current is about $6.5~{}mA$ . In all the experiments reported here
the laser is driven at a constant current of 9 mA. The laser beam is focused
by a high numerical aperture (NA=0.7) lens on the central region of a basic
computer speaker located at about 20 cm from the laser. This speaker is put in
motion via an electrical signal produced by the sound card of a computer which
can easily produce many kinds of patterns at a sampling rate of 44.1 kHz. The
linearity of the speaker response has been assessed over a range of
frequencies from 5 Hz to 100 Hz and over the range of voltage provided by the
sound card. Over this range, the speaker responds with constant 0 phase. Thus,
in the range where the linearity of the displacement has been checked, one can
use the voltage at the speaker as an independent measurement of the position
of the target with respect to some unknown origin. From that point on, we will
therefore use this voltage as a proxy for the target position. The self mixing
signal is measured as a voltage at the laser electrodes, which is amplified by
an AC-coupled amplifier with $10^{4}$ amplification factor and several MHz
bandwidth. We deliberately did not optimize the self-mixing signal quality,
exactly because one of our aims is to check that neural networks can help in
making the measurement work even in sub optimal conditions.

###  2.2 Network setup

The first thing to be noted is that in the "few wavelengths" range of
displacement, the self-mixing signal contains no information about the
absolute position. Although one can be tempted to consider that counting
fringes (or some other equivalent technique) will lead to knowledge of the
exact position with respect to some unknown arbitrary origin, this approach is
bound to diffuse in the long term: If for some reason a fringe is missed, the
measurement system has no way to recover from this error because the physics
of the system does not include this information. Thus, (independently of how
accurate the fringe counting is unless it is strictly perfect), a position
measurement will unavoidably loose accuracy at a rate proportional to
$\sqrt{t}$ where $t$ is the measurement duration. Therefore, our aim here is
to provide a measurement of the displacement within a prescribed time
interval, a velocity.

Once this is established, the setting of the architecture of the neural
network is strongly influenced by the specific question one wants to address.
Here we assume that the self-mixing signal is acquired at a much larger
sampling rate than the Nyquist frequency of the displacement signal one wants
to measure. This is a very reasonable requirement in this context since most
approaches address the question by counting fringes. Here we assume that the
signal is sampled at least 256 times faster than the Nyquist frequency of the
displacement to be measured. Therefore, the reconstruction of the trajectory
consists in inferring from 256 self-mixing signal points one single
instantaneous velocity corresponding to the displacement of the target during
the 256 points acquisition. Then, in terms of machine learning the problem is
reduced to a "regression" problem, where some algorithm must provide a single
number on the basis of the available information (a piece of time trace of
length 256).

The setting therefore consists in analyzing a sequence where temporal ordering
matters and therefore a recurrent neural network can be envisioned as a
suitable architecture. However, these networks are notoriously difficult to
train and convolutional neural networks are known to be an easier to train and
valid alternative alternative. Therefore, we build a network based on a stack
of 1-dimensional convolutional layers with pooling layers between two
convolutional layers. At the end of the stack, two fully connected layers
convert the features identified by the convolutional layers into a single
number which is the inferred velocity of the target during the measurement
sequence. More details are given in appendix, table 1. This global
architecture was chosen from first principles of neural network design [15]
and the model details where then determined empirically. The network was
implemented with the Keras library, which offers an excellent tradeoff in
terms of complexity and versatility for our purpose [22].

###  2.3 Network training

Once the network architecture is chosen, the network must be trained with
known data. In practice, that means providing the network a large number of
pairs $[s(t_{0},...,t_{0}+256~{}dt),v]$ where $s(t)$ is a self-mixing signal
acquired during 256 sampling times $dt$ and $v$ is the average velocity of the
target during the duration of the interval $256~{}dt$. One must underline that
neural networks are known to be able to represent arbitrary functions provided
a sufficient number of layers and cells are present in the network [23].
Therefore, given enough computer time for training, a sufficiently large
network will be able to perfectly reproduce the training data it has been
shown. This means that a model trained this way achieves excellent accuracy.
However, one of the key issues with self-mixing implementations is that the
alignment conditions are sometimes different from one measurement to the next.
Equivalently, speckle generated by the reflection of the beam on a diffusive
target will lead to effective variations of the feedback strength parameter in
the course of the measurement. Therefore, the network trained here must be
able to adapt to these changes. This is known as the capacity of the network
to generalize the features which were learnt and identify them in unseen data.
To achieve this, we train the network on a deliberately limited set of data
and observe the reconstruction of the network on a very different data set,
both in terms of the dynamics of the target (different displacement patterns)
and in terms of the alignment of the beam on the target. The training data
consists exclusively of measurements of the interferometric signal in response
to periodic displacement of the target. We record self mixing signals in six
different alignment conditions, in three of them a double peak is visible in
each fringe. For each of these alignment conditions, we record the self mixing
signal for a set of 19 frequencies evenly spaced between 10 and 100Hz. For
each frequency we record 5 different amplitude signals. For each of these
settings we record sinusoidal and triangular waveforms. The sampling rate of
the oscilloscope is set to 250 kHz so that $dt=4\mu$s. In total the network is
trained on about $1.95\times 10^{3}$ segments of 256 time steps, each of them
of duration $256*4\mu$s $=1.024$ ms. The operation regime of self-mixing
sensing setups is often characterized in terms of the $C$ feeback parameter.
Here, the alignment configurations we use are such that the system operates in
the weak feedback regime $C<1$ (we do not observe multistability). However, we
also avoid the weakest feeback regime $C<<1$ in which the interferometric
signal is symmetric since it does not carry the relevant information. As is
common in deep learning network training, the data is further augmented by
adding noise to the training set. Here we add a delta-correlated gaussian
noise on top of the measured interferometric signal. It is trained by
minimizing the mean squared error between a guess it provides and the known
measured displacement.

##  3 Results

After training, one will assess the performance of the neural network (also
"the model") by comparing the displacement reconstructed from the
interferometric signal and the voltage at the speaker’s ends, used as a proxy
of position. First, we check the model’s prediction accuracy in known settings
(periodic signals and known alignment conditions) and then in unseen settings.

Figure 2: Reconstruction of periodic displacements similar to those used to
train the network. Top row: interferometric signal corresponding to three
different alignment configurations. Bottom row: measured displacement signal
(blue continuous line) and trajectory predicted by the network. The orange
crosses are predicted by the model and the dashed line is a simple cubic
interpolation. The shaded areas correspond to the interferometric signals
shown on the top row.

###  3.1 Periodic signals, known alignment conditions

We show on Fig. 2 how the network reconstructs examples of periodic traces
after training. This exact sequence has not been used during training but
these alignment conditions were used during training and sequences with
identical frequencies and amplitude were used during training.

On the top row we show three examples of interferometric signal which
correspond to three of the alignment conditions used during the training of
the network and three different displacement frequencies. On the bottom row,
we show the displacement per time unit of the target, as it can be measured
from the voltage at the edges of the speaker (blue continuous line).
Independently of that voltage measurement, we use the trained neural network
to infer the displacement from the self-mixing signal. This is the orange
dashed line, which is almost perfectly superimposed to the actual displacement
measured from the voltage at the speaker’s ends. This almost perfect
reconstruction is not very surprising since, even if the network had not seen
this exact piece of time trace during training, it has seen periodic signals
at these frequencies, these amplitudes and in these exact alignment
conditions. That is however a confirmation that the training of the network
has worked to an excellent accuracy and under different alignment conditions.

###  3.2 Aperiodic signals, unknown alignments

We check the capacity of the statistical model to adapt to unseen situations
by preparing a completely different displacement pattern. This pattern is
obtained by applying a fifth order butterworth band-pass filter between 5 and
100 Hz to a delta-correlated gaussian random noise. This pattern is sent to
the speaker in two different alignment conditions, none of them corresponding
the the situations used during training. In one of the two situations, the
interferometric signal shows a double-peak structure. The two interferometric
signals are then concatenated into a single time series and we use the model
to reconstruct the displacement of the target corresponding to this
concatenated time series. The results are shown on Fig. 3.

Figure 3: Reconstruction of unknown and complex displacements in alignment
situations which are not in the training set. Top: the blue line is the
displacement measured from the speaker voltage, the orange line is the model’s
prediction. Middle row: zoom around the central area of the top row. The
discontinuity close to 100 ms is where we numerically connected the two time
traces (see text). Bottom row: interferometric signal corresponding to the
middle row.

As can be immediately appreciated, the reconstruction is excellent, the
prediction matching almost perfectly the independently measured displacement.
Of course the discontinuity close to 100 ms, where the two measurements are
artificially concatenated, cannot be predicted by the network since it is
absent in the measured interferometric signal. We just choose to emphasize
this region as it shows that the prediction is essentially insensitive to the
alignment conditions, which change abruptly in the middle of the trace. As is
evident from the lower panel of Fig. 3, reconstructing a trajectory from this
interferometric signal would be very difficult due to the presence of noise
and very widely varying fringe shapes and repetition rates.

From the above, one concludes that the model is able to generalize from its
learning set to provide an accurate reconstruction of the displacement in
unseen alignment conditions and for very complex time series, much more
difficult to analyze than the simple periodic time traces used during the
training phase.

To better appreciate the accuracy of the inference, one can plot the predicted
displacement as a function of the actual displacement as shown on Fig. 4. A
perfect reconstruction would be the one shown by the orange line where the
prediction is exactly equal to the truth. We can quantify the reconstruction
quality by the Pearson’s correlation coefficient between the reconstruction
and the ground truth which is here 0.90 and the absolute standard error which
is here $0.30\lambda/ms$. Specifically, one can notice that the prediction is
less good for the largest absolute values of displacement. This can be related
to the statistical properties of the training set as shown on the right panel
of Fig. 4. Here one can appreciate that absolute values of displacement larger
than 2.5$\lambda/ms$ have been seen by the network during training only a few
hundreds of times, while smaller displacements are much more frequent in our
training set. Thus, the large displacements are very under-represented in the
training set. This results in a lower precision of the reconstruction for
larger displacements, which can also be appreciated on the top panel of Fig. 3
where the largest displacements are in general under estimated.

Figure 4: The accuracy of prediction on unknown samples (at a given sampling
rate) is conditioned by the training set. Left: displacement inferred by the
model as a function of the true displacement. The correlation is excellent but
the model slightly underestimates the larger displacements (about 2.5
$\lambda$/ms). Right: histograms of the distribution of displacements in the
training sets. The larger displacements (about 2.5 $\lambda$/ms) are very
under-represented in the training set.

###  3.3 Noise sensitivity

One of the difficulties in reconstructing the displacement from the self
mixing signal also comes from the fact that simple Fourier filtering is often
not very efficient at separating the detection noise from the interferometric
signal (although neural networks have been proposed to alleviate this issue
[20]). Here we check that the statistical model is very robust to the addition
of noise on top of the interferometric signal. To assess this robustness, we
use the model to reconstruct the displacement corresponding to the complex
interferometric signal described in 3.2 after adding to this signal a gaussian
white noise. As we show on Fig. 5, the model predictions are extremely robust.
Since the added noise is $\delta-$ correlated, its standard deviation
$\sigma_{n}$ is a measure of its power density. Here, we normalize the
interferometric signal itself in absence of added noise to its standard
deviation $\sigma_{s}$ so that $\sigma_{s}=1$. We then vary the standard
deviation of the added noise $\sigma_{n}$ between 0 and 3 times $\sigma_{s}$.
On Fig. 5a), we observe that both the root mean squared and the absolute error
remain very low up to $\sigma_{n}=1$ where it grows significantly. On Fig.
5b), we plot the correlation coefficient between the reconstructed
displacement and the independently measured displacement. As for the error,
the correlation coefficient indicates an excellent reconstruction of the
displacement up to approximately $\sigma_{n}=1$. We show the prediction and
the noisy interferometric signal for $\sigma_{n}=0.7$ on Fig. 5 c) and d)
respectively. As can be easily observed, the interferometric signal would be
rather difficult to process by the usual means and the stochastic model
provides a very useful reconstruction.

Figure 5: Robustness against detection noise. a) Mean absolute (orange
continuous) and RMS (blue dashed) error of the reconstruction as function of
added noise in units of the interferometric signal’s standard deviation
$\sigma_{n}$. b) Pearson’s correlation coefficient between the prediction and
the signal as function of added noise. c) Example of prediction (orange line)
and true displacement (blue line) for an added noise power density
$\sigma_{n}=0.7\times\sigma_{n}$ where $\sigma_{n}$ is the interferometric
signal’s standard deviation. d) Interferometric signal in the same situation,
zoom over the 94-106 ms region (the same signal as Fig. 3, with added noise).

###  3.4 Unknown experiment

The analysis above has shown that the statistical model is able to reconstruct
the displacement from the self-mixing signal in a broad range of unknown
conditions. However, all of the above was realized on a single experimental
setup. Contrary to a physical model, which is constructed to capture only the
universal features of an experiment, an empirically constructed statistical
model such as the neural network we use may capture also non-universal and
system-specific features. Thus, it is interesting to check what the model can
predict on the basis of a different experiment, based on the same principle.
To address this question, we prepare an "almost-twin" experiment, based on the
same self mixing interferometry principle shown in Fig. 1 but featuring a
different laser (HL6323MG, $\lambda=639~{}nm$, driven at $I=75~{}mA$ for a
threshold current $I_{th}=45~{}mA$), a different speaker (with a different
range of linear response), a different voltage amplifier for the acquisition
of the laser diode voltage etc. Although this experiment is in principle the
same, it differs in many of the details which should not be relevant to the
physics, yet carry a significant risk of distortion of the interferometric
signal as compared to the one used in training set.

To check the model’s ability to analyze this new self mixing experiment, we
prepare a new displacement time series consisting of the sum of four different
frequencies $S=\sum_{i=1}^{4}\sin{2\pi f_{i}t}$ with
$f_{i}=425,718,808,1076~{}Hz$. This time series is therefore in a very
different (much higher) frequency band with respect to the training
experiment. In order to provide the model with comparable input data, the
acquisition of the self-mixing signal is performed at a ten times faster rate
than in the training experiment (2.5 MHz). The displacement per time unit of
the speaker is, as in the previous experiment, measured as a voltage at the
edges of the speaker. The comparison between the displacement estimated from
the speaker’s voltage and the displacement reconstructed from the
interferometric signal is shown on Fig. 6.

Figure 6: Reconstruction of the displacement on an unknown experiment. Top:
comparison between the displacement predicted by the model (orange line) and
the actual displacement (blue line). There are no free parameters. Middle:
zoom on one specific region. Bottom: self mixing signal corresponding to
middle panel. The green circles show where jumps between bistable states
occur.

The agreement between the prediction and the measurement is strikingly good,
especially taking into account that no free parameter exist: The model trained
on experiment 1 can immediately be used to infer displacements in units of
$\lambda/dt$ in experiment 2.

It is important to underline once more the robustness of this process with
respect to specific experimental conditions. For instance, in this experiment,
the interferometric signal shows clear signs of bistability between external
cavity modes in forms of very fast jumps between states (green circles on
bottom panel of Fig. 6). These features are absent from the training set. Here
one sees that the model essentially filters them out automatically. Besides
this, it is also worth noting that, as compared to Fig. 3 for instance, the
displacements and time scales are very different. This shows that, provided an
adequate sampling rate is chosen, the model can work at much higher frequency
than the band it was trained in and for much larger displacements per time
unit. This feature is not unexpected since the network knows only about
displacement "per $256\times dt$", without reference to the exact value of
$dt$. Therefore, with adequate sampling, the measurement range can be
tremendously extended with respect to the training range. This feature is
extremely useful since it allows training in an easily accessible range
(displacement and frequency band) and prediction in very different range for
more demanding applications.

##  4 Discussion

The results above clearly show that a convolutional neural network is a useful
tool in the reconstruction of a target’s displacement, very robust to unknown
displacement signal shapes, alignment conditions, electronic noise and even
whole setups. We have also verified that a network trained in the 10-100 Hz
frequency band can also meaningfully reconstruct displacements including
frequencies of hundreds of Hz, provided the measurement sampling rate is
adapted. Thus, the neural network strength lies less in the absolute precision
that it allows than in its robustness against detailed experimental conditions
and versatility across the "sub-wavelength/analog" and "beyond
wavelength/digital" classification [6] for arbitrarily complex waveforms.

One natural question which arises when preparing a neural network is that of
model capacity [15]. A network which does not possess enough cells or layers
may be unable to take into account all the complexity of the task. On the
other hand, a network with a very large number of cells and layers will sooner
or later learn features of the experiment which should not be significant (for
instance, all the details about an amplifier used in the setup). This prevents
the network from generalizing, ie accurately predicting unknown data. This is
in principle dealt with during the training phase [15] but it is only when the
network processes fully new data that this issue can be totally ruled out.
Here this issue has been taken care of by predicting arbitrarily complex
trajectories and also by using two different setups. In fact, the imperfect
reconstruction in the case of the unknown experiment is most probably due to
the model learning some system-specific features of the training experiment.
This can be mitigated by a minor retraining of the final layer of the model on
the new experiment (a procedure known as "fine tuning" in the deep learning
context). We have noticed that a larger network featuring more than $10^{5}$
coefficients instead of the $5.7\times 10^{4}$ used here does not lead to
better training and may even lead to worse predictions in the unknown
experiment.

The performance of the network is of course strongly related to the training
data set which is used. Here we deliberately use only a very limited set of
displacements during training in order to very clearly show the generalization
phenomenon, the network being able to predict correctly displacement shapes it
has never seen before. For real use beyond the proof of concept presented
here, more refined training is possible: A training set featuring a more
uniform distribution of displacements will provide a more accurate
reconstruction of the larger displacements for a given sampling rate. As an
alternative, simply increasing the sampling rate at the prediction time may
also be a sufficient solution to adapt the time series to the operating range
of the model as we have shown in 3.4. Care must be taken when training the
network that the correct operating range of the neural network is set by a
displacement per time unit, which includes limitations in terms of
displacement frequency and amplitude. Translating it in terms of counting
fringes, that means that the network will saturate beyond a certain number of
fringes during the $256\times dt$ measurement window. In terms of feedback
range, here we have used only $C<1$ in training, avoiding too low values of
$C$ where the interferometric signal is symmetric. At the prediction phase,
the model is robust to $C$ slightly overcoming unity but when multistability
becomes strong the information of few wavelengths displacements is lost and
the model has no chance to recover it. Similarly, we have checked that when
$C$ is so low that the signal is symmetric, the model cannot predict
accurately. A full characterization of performance degradation and the use of
multichannel measurements to mitigate this issue is beyond the scope of this
work.

One particularly interesting avenue to circumvent the limitations of an
experimental training set is to train the network on numerically generated
data [24]. One drawback is that the network will not learn more than what is
in the physical model used for the simulations, which may be hard in complex
settings such as multimode lasers [25, 26]. However, numerics may provide a
way to obtain a training set for which controlled laboratory experiments would
be very hard to realize such as hard shocks or high frequency and high
amplitude displacements. Experimental data may then be used to refine the
training by using different sampling times as described in the previous
paragraph.

Once trained, a convolutional neural network can be used in real time since no
pre-processing of the data is required and prediction over thousands of
interferometric measurements is very fast: As an example, 3 seconds of signal
(749056 interferometric data points) are processed in 0.16 seconds on a
standard laptop. In addition, after the initial training, a neural network can
relatively easily be repurposed by retraining only its final layers even with
a very limited set of data. For instance, it would be particularly interesting
to assess the performance of the network trained here on self-mixing schemes
which include bias current modulation towards some other sensing task such as
refractive index measurements. Alternatively, the input layer of the network
can also be reworked at minor cost to take into account multichannel
measurements and most interesting would probably be to integrate this approach
into multimodality imaging systems [27].

##  5 Conclusion

To conclude, we have presented a detailed analysis of how a reasonably simple
convolutional neural network can be used to reconstruct the displacement of a
target on the basis of self mixing interferometry. We believe that this
approach can become one of the many tools which can be used to tailor or
enhance self-mixing coherent sensing setups. Far from being limited to
displacement measurement and single mode settings, we believe that computer
neural networks can become an extremely useful element of many sensing
apparatus, especially self-mixing setups. Finally, we stress that there are
very few hard rules about the design of neural networks and this design is in
itself often an area of research. The architecture we use here is essentially
a simple starting point and many refinements are possible. More specifically,
one of the most immediate extensions of the work exposed above is to train a
network which can provide an estimation of the accuracy of the reconstruction.
This can be achieved by adding a calibrated regression stage on top of the
convolutional base prepared here [28]. Other extensions may include more
complex network topologies, perhaps mixing convolutional and recurrent layers
or including skip connections.

##  Appendix A Network details

The key elements of the network are shown on table 1. We refer the reader to
deep learning fundamentals for background information [15]. The total number
of parameters (57 153) is much smaller than the number of time series segments
in the data set even before augmentation. Networks of identical architecture
with more cells per layer did not lead to significant improvements. A dropout
layer is used here mostly as "safety net" since the data set is very large
anyway which makes overfitting improbable. The training of the network takes
about ten minutes on a simple GPU (GeForce GTX 1060) and about four times more
on CPU (Intel Xeon 3.8GHz).

Network structure: sequential
---
Layer type |  Main
hyperparameters |  Trainable
parameters |  Output shape
1D convolutional |  kernel size: 7
filters: 16 |  128 |  (250, 16)
Max Pooling |  pool size: 2 |  0 |  (125, 16)
1D convolutional |  kernel size:7
filters: 32 |  3616 |  (119, 32)
Max Pooling |  pool size: 2 |  0 |  (59, 32)
1D convolutional |  kernel size: 7
filters 64 |  14400 |  (53, 64)
Max Pooling |  pool size: 2 |  0 |  (26, 64)
Dropout |  10% |  0 |  (26, 64)
1D convolutional |  kernel size: 7
filters: 64 |  28736 |  (20, 64)
Max Pooling |  pool size: 2 |  0 |  (10, 64)
Fully connected |  units: 16 |  10256 |  16
Fully connected |  units: 1 |  17 |  1
Table 1: Main parameters of the network used in this work. The network is a
sequence of 1-dimensional convolutional and dropout layers followed by two
fully connected layers for the final regression. The total number of trainable
parameters is 57 153.

## Acknowledgments

The authors thank Dr. L. Columbo, Dr. M. Dabbicco and Dr. F. Pedaci for
helpful discussions.

## Disclosures

The authors declare no conflicts of interest.

## References

  * [1] G. Giuliani, M. Norgia, S. Donati, and T. Bosch, “Laser diode self-mixing technique for sensing applications,” Journal of optics A: Pure and applied optics 4, S283 (2002).
  * [2] D. M. Kane and K. A. Shore, _Unlocking dynamical diversity: optical feedback effects on semiconductor lasers_ (John Wiley & Sons, 2005).
  * [3] S. Donati, “Developing self-mixing interferometry for instrumentation and measurements,” Laser & Photonics Reviews 6, 393–417 (2012).
  * [4] T. Taimre, M. Nikolić, K. Bertling, Y. L. Lim, T. Bosch, and A. D. Rakić, “Laser feedback interferometry: a tutorial on the self-mixing effect for coherent sensing,” Advances in Optics and Photonics 7, 570–631 (2015).
  * [5] J. Li, H. Niu, and Y. X. Niu, “Laser feedback interferometry and applications: a review,” Optical Engineering 56, 050901 (2017).
  * [6] S. Donati and M. Norgia, “Overview of self-mixing interferometer applications to mechanical engineering,” Optical Engineering 57, 051506 (2018).
  * [7] R. C. Addy, A. W. Palmer, K. Thomas, and V. Grattan, “Effects of external reflector alignment in sensing applications of optical feedback in laser diodes,” Journal of lightwave technology 14, 2672–2676 (1996).
  * [8] L. Lv, H. Gui, J. Xie, T. Zhao, X. Chen, A. Wang, F. Li, D. He, J. Xu, and H. Ming, “Effect of external cavity length on self-mixing signals in a multilongitudinal-mode fabry–perot laser diode,” Applied optics 44, 568–571 (2005).
  * [9] M. Norgia, S. Donati, and D. D’Alessandro, “Interferometric measurements of displacement on a diffusing target by a speckle tracking technique,” IEEE Journal of quantum electronics 37, 800–806 (2001).
  * [10] U. Zabit, R. Atashkhooei, T. Bosch, S. Royo, F. Bony, and A. Rakic, “Adaptive self-mixing vibrometer based on a liquid lens,” Optics letters 35, 1278–1280 (2010).
  * [11] O. D. Bernal, U. Zabit, and T. M. Bosch, “Robust method of stabilization of optical feedback regime by using adaptive optics for a self-mixing micro-interferometer laser displacement sensor,” IEEE Journal of Selected Topics in Quantum Electronics 21, 336–343 (2014).
  * [12] A. L. Arriaga, F. Bony, and T. Bosch, “Speckle-insensitive fringe detection method based on hilbert transform for self-mixing interferometry,” Applied optics 53, 6954–6962 (2014).
  * [13] M. Usman, U. Zabit, O. D. Bernal, G. Raja, and T. Bosch, “Detection of multimodal fringes for self-mixing-based vibration measurement,” IEEE Transactions on Instrumentation and Measurement 69, 258–267 (2019).
  * [14] A. A. Siddiqui, U. Zabit, O. D. Bernal, G. Raja, and T. Bosch, “All Analog Processing of Speckle Affected Self-Mixing Interferometric Signals,” IEEE Sensors Journal 17, 5892–5899 (2017).
  * [15] I. Goodfellow, Y. Bengio, A. Courville, and Y. Bengio, _Deep learning_ , vol. 1 (MIT press Cambridge, 2016).
  * [16] H. Li, C. Zhang, N. Song, and H. Li, “Deep learning-based interference fringes detection using convolutional neural network,” IEEE Photonics Journal 11, 1–14 (2019).
  * [17] A. Reyes-Figueroa and M. Rivera, “Deep neural network for fringe pattern filtering and normalisation,” arXiv preprint arXiv:1906.06224 (2019).
  * [18] K. Kou, C. Wang, T. Lian, and J. Weng, “Fringe slope discrimination in laser self-mixing interferometry using artificial neural network,” Optics & Laser Technology 132, 106499 (2020).
  * [19] L. Wei, J. Chicharo, Y. Yu, and J. Xi, “Pre-processing of signals observed from laser diode self-mixing intereferometries using neural networks,” in _2007 IEEE International Symposium on Intelligent Signal Processing,_ (IEEE, 2007), pp. 1–5.
  * [20] I. Ahmed, U. Zabit, and A. Salman, “Self-mixing interferometric signal enhancement using generative adversarial network for laser metric sensing applications,” IEEE Access 7, 174641–174650 (2019).
  * [21] X.-l. Wang, L.-p. Lü, L. Hu, and W.-c. Huang, “Real-time human blood pressure measurement based on laser self-mixing interferometry with extreme learning machine,” Optoelectronics Letters 16, 467–470 (2020).
  * [22] F. Chollet _et al._ , “keras,” (2015).
  * [23] K. Hornik, M. Stinchcombe, H. White _et al._ , “Multilayer feedforward networks are universal approximators.” Neural networks 2, 359–366 (1989).
  * [24] R. Kliese, T. Taimre, A. A. A. Bakar, Y. L. Lim, K. Bertling, M. Nikolić, J. Perchoux, T. Bosch, and A. D. Rakić, “Solving self-mixing equations for arbitrary feedback levels: a concise algorithm,” Applied optics 53, 3723–3736 (2014).
  * [25] L. Columbo, M. Brambilla, M. Dabbicco, and G. Scamarcio, “Self-mixing in multi-transverse mode semiconductor lasers: model and potential application to multi-parametric sensing,” Optics express 20, 6286–6305 (2012).
  * [26] L. Columbo and M. Brambilla, “Multimode regimes in quantum cascade lasers with optical feedback,” Optics Express 22, 10105–10118 (2014).
  * [27] M. Brambilla, L. L. Columbo, M. Dabbicco, F. De Lucia, F. P. Mezzapesa, and G. Scamarcio, “Versatile multimodality imaging system based on detectorless and scanless optical feedback interferometry—a retrospective overview for a prospective vision,” Sensors 20, 5930 (2020).
  * [28] V. Kuleshov, N. Fenner, and S. Ermon, “Accurate uncertainties for deep learning using calibrated regression,” in _Proceedings of the 35th International Conference on Machine Learning,_ vol. 80 of _Proceedings of Machine Learning Research_ J. Dy and A. Krause, eds. (PMLR, Stockholmsmässan, Stockholm Sweden, 2018), pp. 2796–2804.