Laurent Bétermin

Ladislav Šamaj

Igor Travěnec

On off-critical zeros of lattice energies in the neighborhood of the Riemann zeta function

Keywords: Riemann zeta function, Hurwitz zeta function, critical and off-critical zeros, Riemann hypothesis

The Riemann zeta function ζ(s) := ∞ n=1 1/n s can be interpreted as the energy per point of the lattice Z, interacting pairwisely via the Riesz potential 1/r s . Given a parameter ∆ ∈ (0, 1], this physical model is generalized by considering the energy per point E(s, ∆) of a periodic one-dimensional lattice alternating the distances between the nearest-neighbour particles as 2/(1 + ∆) and 2∆/(1 + ∆), keeping the lattice density equal to one independently of ∆. This energy trivially satisfies E(s, 1) = ζ(s) at ∆ = 1, it can be easily expressed as a combination of the Riemann and Hurwitz zeta functions, and extended analytically to the punctured s-plane C \ {1}. In this paper, we perform numerical investigations of the zeros of the energy {ρ = ρ x + iρ y }, which are defined by E(ρ, ∆) = 0. The numerical results reveal that in the Riemann limit ∆ → 1 -theses zeros include the anticipated critical zeros of the Riemann zeta function with ℜ(ρ x ) = 1 2 as well as an unexpected -comparing to the Riemann Hypothesis -infinite series of off-critical zeros. The analytic treatment of these off-critical zeros shows that their imaginary components are equidistant and their real components diverge logarithmically to -∞ as ∆ → 1 -, i.e., they become invisible at the Riemann's ∆ = 1.

Introduction and main results

Let two points at distance r interact via the Riesz potential 1/r s with real s [START_REF] Brauchart | Optimal discrete Riesz energy and discrepancy[END_REF]. If the points are located on the lattice Z and interact pairwisely by the Riesz potential where s > 1, the energy per point is given by the Riemann zeta function [START_REF] Riemann | Über die Anzahl der Primzahlen unter einer gegebenen Grösse[END_REF] ζ(s

) := 1 2 n∈Z * 1 |n| s = ∞ n=1 1 n s s > 1, (1)
where the prefactor 1 2 is due to the fact that each interaction energy is shared by a pair of points. The function ζ can be analytically continued to the whole complex s-plane, with a simple pole at s = 1. The Riemann zeta function plays a fundamental role in the algebraic and analytic number theories [START_REF] Hadamard | Étude sur les propriétés des fonction entières et un particulier d'une fonction considéré par Riemann[END_REF][START_REF] Hardy | Sur les zeros de la fonction ζ(s)[END_REF][START_REF] Riesz | Sur l'hypothèse de Riemann[END_REF][START_REF] Hardy | The zeros of Riemann's zeta-function on the critical line[END_REF][START_REF] Hutchinson | On the Roots of the Riemann Zeta-Function[END_REF][START_REF] Titchmarsh | The Zeros of the Riemann Zeta-Function[END_REF][START_REF] Selberg | Contributions to the theory of the Riemann zeta-function[END_REF], see monographs [START_REF] Edwards | Riemann's Zeta function[END_REF][START_REF] Ivić | The Riemann Zeta Function[END_REF][START_REF] Titchmarsh | The Theory of The Riemann Zeta-function[END_REF]. The so-called Riemann Hypothesis about the location of its nontrivial zeros exclusively on the critical line ℜ(s) = 1 2 (the symbol ℜ means the real part) is one of the Hilbert and Clay Millennium Prize problems [START_REF] Jaffe | The Millenium Grand Challenge in Mathematics[END_REF]. Throughout the present paper we assume that the Riemann Hypothesis holds. The Riemann zeta function and its Epstein's [START_REF] Epstein | Zur Theorie allgemeiner Zetafunctionen[END_REF][START_REF] Epstein | Zur Theorie allgemeiner Zetafunctionen II[END_REF][START_REF] Chowla | On Epstein's zeta function[END_REF][START_REF] Travěnec | Generation of off-critical zeros hypercubic Epstein zeta functions[END_REF]], Hurwitz's [START_REF] Hurwitz | Einige Eigenschaften der Dirichletschen Fuctionen F (s) = D n • 1 n s die bei der Bestimmung der Klassenzahlen binärer quadratischer Formen auftreten[END_REF][START_REF] Spira | Zeros of Hurwitz zeta functions[END_REF][START_REF] Nakamura | Real zeros of Hurwitz-Lerch zeta and Hurwitz-Lerch type of Euler-Zagierdouble zeta functions[END_REF]], Barnes's [START_REF] Barnes | On the theory of the multiple gamma function[END_REF], etc. generalisations have numerous applications both in mathematics (prime numbers, applied statistics [START_REF] Borwein | Lattice sums then and now[END_REF]) and in physics [START_REF] Elizalde | Ten Physical Applications of Spectral Zeta Functions[END_REF].

Let the Riemann zeta function be a member of a family of functions which exhibit nontrivial zeros off the critical line. Possible mechanisms of the disappearance of these off-critical zeros at the Riemann's point might of general interest. In this paper, we propose a natural extension of the Riemann zeta function as the energy of a unit density lattice L ∆ with alternating distances between the nearest neighbours, say 2/(1 + ∆) and 2∆/(1 + ∆); due to the ∆ → 1/∆ symmetry of the problem, it is sufficient to restrict oneself to ∆ from the interval (0, 1].

In analogy with the original model with constant unit spacing, each point interacts pairwisely with the other points via the Riesz interaction 1/r s , s > 1 and the lattice energy per point is therefore given (see Proposition 2.1) by

E(s, ∆) = 1 2 s ζ(s) + 1 2 s+1 ζ s, 1 1 + ∆ + ζ s, ∆ 1 + ∆ , s > 1, ∆ ∈ (0, 1],
where

ζ(s, a) = ∞ n=0 1 (n + a) s , s > 1, (2)
is the Hurwitz zeta function with the real (positive) parameter a. Remark that this lattice energy, as a combination of Riemann and Hurwitz zeta functions, has an analytic continuation on C\{1} (see Proposition 2.3).

For given ∆ ∈ (0, 1], the set of zeros of the lattice energy is defined as

Z ∆ := {ρ = ρ x + iρ y ∈ C, (ρ x , ρ y ) ∈ R 2 : E(ρ, ∆) = 0}, noticing that Z 1
is the set of zeros of the Riemann zeta function composed by trivial zeros (i.e., ρ ∈ -2N) and critical zeros (i.e., ℜ(ρ) = 1 2 ) assuming that Riemann Hypothesis holds. Furthermore, for specific values of the parameter ∆ ∈ {1/5, 1/3, 1/2}, the energy can be factorized as

E(s, ∆) = f ∆ (s)ζ(s)
where f ∆ is a sum of p s with integers p. This automatically gives us critical zeros (i.e. solutions of ζ(ρ) = 0 assuming the Riemann Hypothesis) and possible off-critical zeros (i.e. solutions of f ∆ (ρ) = 0), as shown in Proposition 2.4.

The goal of this paper is to study, both numerically and analytically, the set of zeros Z ∆ when ∆ is in a neighborhood of 1, i.e. when E(•, ∆) is in the neighborhood of the Riemann zeta function. Numerical and analytic analysis shows that approaching ∆ → 1 -the zeros of E(ρ, ∆) involve the anticipated critical zeros of the Riemann zeta function with ℜ(ρ x ) = 1 2 as well as an infinite series of unexpected off-critical zeros with the following asymptotics for their components, as ∆ → 1 -(see Theorem 4.2):

ρ x (∆) = 2 ln 2 ln(1 -∆) + O(1 -∆) → -∞, ρ y (∆) = (2k + 1)π ln 2 + O (1 -∆) 2 ln 3 ln 2 , k ∈ Z.
This means that, asymptotically, there is a infinite sequence of equidistant zero components along the ρ y axis. Furthermore, the divergence of ρ x to -∞ as ∆ → 1 -is an example of the disappearance of off-critical zeros when approaching the Riemann's point. Moreover, the behavior of these zero components with respect to ∆ ∈ (0, 1] is numerically studied (see Figures 1 and2).

Plan of the paper. The generalized 1D model for Riesz points with alternating lattice spacings is presented in section 2. The energy per particle E(s, ∆) is expressed as a combination of Hurwitz zeta functions in section 2.1. The properties of the Hurwitz zeta function are discussed in section 2.2. Special values of the parameters ∆ when the energy E(s, ∆) factorizes itself onto the product of the Riemann zeta function and some simple function are given in section 2.3. Numerical results for zeros at any 0 < ∆ < 1, together with tests at the special values of ∆ = 1/5, 1/3, 1/2 are presented in section 3. The spectrum of critical and off-critical zeros in the Riemann's limit ∆ → 1 -is discussed in section 4.

The generalized one-dimensional model

Definition of the model

Given ∆ ∈ (0, 1], we consider the infinite set of points L ∆ ⊂ R given by

L ∆ := 2Z ∪ 2Z + 2∆ 1 + ∆ ,
which is the unit density periodic configuration with alternating distances 2/(1 + ∆) and 2∆/(1 + ∆), since 2 -2/(1 + ∆) = 2∆/(1 + ∆). Assuming that each pair of points in L ∆ interacts via the Riesz potential 1/r s , s > 1, the total energy per point of this system is therefore

E(s, ∆) := 1 4 k∈{0, 2∆ 1+∆ } p∈L ∆ p =k 1 |p -k| s .
The following proposition shows how to write this energy in terms of Riemann and Hurwitz zeta functions.

Proposition 2.1. For any s > 1 and any ∆ ∈ (0, 1], we have

E(s, ∆) = 1 2 s ζ(s) + 1 2 s+1 ζ s, 1 1 + ∆ + ζ s, ∆ 1 + ∆ . (3)
Proof. We simply compute the above double sum as follows:

E(s, ∆) : = 1 4 k∈{0, 2∆ 1+∆ } p∈L ∆ p =k 1 |p -k| s = 1 4 k∈{0, 2∆ 1+∆ }     p∈2Z p =k 1 |p -k| s + p∈2Z+ 2∆ 1+∆ p =k 1 |p -k| s     = 1 4 n∈Z n =0 1 |2n| s + 1 4 n∈Z 2n = 2∆ 1+∆ 1 2n -2∆ 1+∆ s + 1 4 n∈Z 2n =-2∆ 1+∆ 1 2n + 2∆ 1+∆ s + 1 4 n∈Z n =0 1 |2n| s = 1 2 s+1 n∈Z * 1 |n| s + 1 2 s+2 n∈Z 1 n + ∆ 1+∆ s + 1 2 s+2 n∈Z 1 n -∆ 1+∆ s .
We now split the two last sums in order to get two Hurwitz zeta functions and two rests that we write again in terms of the same Hurwitz zeta functions:

E(s, ∆) : = 1 2 s ζ(s) + 1 2 s+2 ζ s, 1 1 + ∆ + ζ s, ∆ 1 + ∆ + 1 2 s+2 +∞ n=1 1 -n + ∆ 1+∆ s + +∞ n=1 1 -n + 1 1+∆ s = 1 2 s ζ(s) + 1 2 s+2 ζ s, 1 1 + ∆ + ζ s, ∆ 1 + ∆ + 1 2 s+2 +∞ n=1 1 n -∆ 1+∆ s + +∞ n=1 1 n -1 1+∆ s .

Since we have, by the change of variables

n = k + 1, +∞ n=1 1 n -∆ 1+∆ s = +∞ k=0 1 k + 1 1+∆ s and +∞ n=1 1 n -1 1+∆ s = +∞ k=0 1 k + ∆ 1+∆ s ,
we obtain

E(s, ∆) = 1 2 s ζ(s) + 2 2 s+2 ζ s, 1 1 + ∆ + ζ s, ∆ 1 + ∆ = 1 2 s ζ(s) + 1 2 s+1 ζ s, 1 1 + ∆ + ζ s, ∆ 1 + ∆
and the proof is complete.

Notice that the energy satisfies the required symmetry relation

E(s, ∆) = E(s, 1/∆).
Remark 2.2 (Crystallization result). It is known (see e.g. [START_REF] Ventevogel | On the configuration of systems of interacting particles with minimum potential energy per particle[END_REF]), by a convexity argument (or by the so-called one-dimensional "Universal Optimality" of Z, see [START_REF] Cohn | Universally optimal distribution of points on spheres[END_REF]) that, for all s > 0, min

∆∈(0,1] E(s, ∆) = E(s, 1) = ζ(s),
with equality if and only if ∆ = 1. From our results (see Theorem 4.2), the Riemann zeta function is therefore at the same time the minimal value of our energy and the only one in its ∆-neighborhood for which the nontrivial zeros are strictly located on the critical line ℜ(s) = 1 2 . It might be interesting to investigate other lattice energies to understand how universal this phenomenon is.

The Hurwitz zeta function and the analytic continuation of E(•, ∆)

The Hurwitz zeta function ( 2) is a generalization of the Riemann zetafunction (1) via a shift a > 0. In particular,

ζ(s, 1) = ζ(s). (4)
In the symbolic computer language Mathematica, the Riemann and Hurwitz zeta functions are tabulated under the symbols Zeta[s] and Zeta[s, a], respectively. The Hurwitz zeta function satisfies two easily verifiable important equalities,

∀x ∈ 0, 1 2 , ∀s > 1, ζ(s, x) + ζ(s, 1/2 + x) = 2 s ζ(s, 2x), (5)
and the multiplication theorem

k ∈ N, ∀s > 1, k s ζ(s) = k n=1 ζ s, n k . (6)
The second relation easily implies that

ζ(s, 1/3) + ζ(s, 2/3) = (3 s -1)ζ(s) (7)
as well as

ζ(s, 1/2) = (2 s -1)ζ(s). (8)
From ( 5) with x = 1/4 and ( 8), we therefore obtain

ζ(s, 1/4) + ζ(s, 3/4) = (4 s -2 s )ζ(s). (9)
From ( 3) and ( 8), it is also straightforward to check that E(s, 1) = ζ(s).

Furthermore, it is clear that (3) holds for all s ∈ C such that ℜ(s) > 1 and we get the following result by directly applying the classical one (see e.g. [START_REF] Fine | Note on the Hurwitz zeta-function[END_REF][START_REF] Hurwitz | Einige Eigenschaften der Dirichletschen Fuctionen F (s) = D n • 1 n s die bei der Bestimmung der Klassenzahlen binärer quadratischer Formen auftreten[END_REF][START_REF] Travěnec | Generation of off-critical zeros hypercubic Epstein zeta functions[END_REF]) on the analytic continuation of the Riemann and Hurwitz zeta functions.

Proposition 2.3. For all ∆ ∈ (0, 1], the function s → E(s, ∆) has an analytic continuation on C\{1}. Furthermore, we have, for all s such that ℜ(s) < 1 and all ∆ ∈ (0, 1],

π -s 2 Γ s 2 E(s, ∆) = 1 2 s+1 ∞ 0 ϑ 1 1 + ∆ , it -1 t -1-s 2 dt + f (s) ,
where ϑ(z, it) = n∈Z e -πn 2 t e 2iπnz is the Jacobi theta function defined for t > 0 and z ∈ C, and

f (s) = ∞ 0 ϑ(0, it) -1 - 1 √ t t s 2 -1 dt for 0 < ℜ(s) < 1, = ∞ 0 ϑ(0, it) - 1 √ t t s 2 -1 dt for ℜ(s) < 0.
Proof. Recall that, for ∆ ∈ (0, 1] and ℜ(s) > 1, we have

E(s, ∆) = 1 2 s ζ(s) + 1 2 s+1 ζ s, 1 1 + ∆ + ζ s, ∆ 1 + ∆ .
It has been shown in [START_REF] Fine | Note on the Hurwitz zeta-function[END_REF][START_REF] Hurwitz | Einige Eigenschaften der Dirichletschen Fuctionen F (s) = D n • 1 n s die bei der Bestimmung der Klassenzahlen binärer quadratischer Formen auftreten[END_REF][START_REF] Travěnec | Generation of off-critical zeros hypercubic Epstein zeta functions[END_REF] that s → ζ(s) and s → ζ(s, a), a > 0, admit an analytic continuation to C\{1}, which implies the same for s → E(s, ∆). Furthermore, writing z = 1 1+∆ , we have

E(s, ∆) = 1 2 s ζ(s) + 1 2 s+1 [ζ (s, z) + ζ (s, 1 -z)] .
Considering the analytic continuation of ζ(1 -α, a), a > 0, the following formula is well-known [START_REF] Fine | Note on the Hurwitz zeta-function[END_REF][START_REF] Hurwitz | Einige Eigenschaften der Dirichletschen Fuctionen F (s) = D n • 1 n s die bei der Bestimmung der Klassenzahlen binärer quadratischer Formen auftreten[END_REF] for all α such that ℜ(α) > 0:

π -1-α 2 Γ 1 -α 2 [ζ (1 -α, z) + ζ (1 -α, 1 -z)] = ∞ 0 [ϑ(z, it) -1]t α 2 dt t (10)
holding for z / ∈ Z (see the remark in Theorem 12.6 on page 257 of [START_REF] Apostol | Introduction to analytic number theory[END_REF]). Replacing α by s = 1 -α (so that ℜ(s) < 1 for ℜ(α) > 0), we write

π -s 2 Γ s 2 E(s, ∆) = π -s 2 2 s Γ s 2 ζ(s) + π -s 2 2 s+1 Γ s 2 [ζ (s, z) + ζ (s, 1 -z)] = 2π -s 2 2 s+1 Γ s 2 ζ(s) + π -s 2 2 s+1 Γ s 2 [ζ (s, z) + ζ (s, 1 -z)] . Next we replace π -s 2 Γ s 2 [ζ (s, z) + ζ (s, 1 -z)
] by the integral given in (10) and 2π -s 2 Γ s 2 ζ(s) by the d = 1 integral in Eq. ( 21) of [START_REF] Travěnec | Generation of off-critical zeros hypercubic Epstein zeta functions[END_REF] for 0 < ℜ(s) < 1 and by the d = 1 integral in Eq. ( 22) of [START_REF] Travěnec | Generation of off-critical zeros hypercubic Epstein zeta functions[END_REF] for ℜ(s) < 0 to complete the proof.

Therefore, we can consider the zeros of s → E(s, ∆) in C\{1} defined as

Z ∆ := {ρ = ρ x + iρ y ∈ C, (ρ x , ρ y ) ∈ R 2 : E(ρ, ∆) = 0},
noticing that Z 1 is the set of zeros of the Riemann zeta function. We recall that, according to the Riemann Hypothesis,

Z 1 = -2N ∪ Z C , Z C ⊂ {ℜ(z) = 1/2},
where Z C is called the set of critical zeros of ζ and -2N are the trivial zeros of ζ. We are going to see in the next sections that Z ∆ can have other nontrivial off-critical zeros.

Factorization and zeros of the energy for special values of ∆

There exist special values of ∆ for which the energy E(s, ∆) factorizes itself into a product of the Riemann zeta function ζ(s) and some simple functions of s, by using the previously presented relations ( 5) and [START_REF] Edwards | Riemann's Zeta function[END_REF]. For these cases, both critical and off-critical zeros can be found easily. The most obvious choice of ∆ is ∆ = 1 for which we have

E(s, 1) = ζ(s). (11)
In the cases ∆ ∈ 1 2 , 1 3

we have the following result giving the zeros of E(s, ∆) as well as the factorization of the energy.

Proposition 2.4. For all s ∈ C\{1}, we have

E(s, 1/2) = 1 2 s+1 (1 + 3 s )ζ(s) and E(s, 1/3) = 1 2 s+1 (2 -2 s + 4 s )ζ(s). (12)
Furthermore, the zeros of E(s, 1/2) and E(s, 1/3) are

Z 1 2 = Z 1 ∪ (2k + 1)iπ ln 3 k∈Z , Z 1 3 = Z 1 ∪ 1 ln 2 ln 1 ± i √ 7 2 + 2πik k∈Z . Proof. For ∆ = 1 2 , one has 1 1 + ∆ = 2 3 and ∆ 1 + ∆ = 1 3
and therefore, using [START_REF] Edwards | Riemann's Zeta function[END_REF], it holds that

E(s, 1/2) = 1 2 s+1 (1 + 3 s )ζ(s). ( 13
)
where the function 1 + 3 s yields an infinite sequence of (purely imaginary) off-critical zeros

ρ k = (2k + 1)iπ ln 3 , k ∈ Z. (14)
Furthermore, for ∆ = 1 3 , we have

1 1 + ∆ = 3 4 and ∆ 1 + ∆ = 1 4
and therefore, applying [START_REF] Epstein | Zur Theorie allgemeiner Zetafunctionen[END_REF], it holds that

E(s, 1/3) = 1 2 s+1 (2 -2 s + 4 s )ζ(s). ( 15
)
The function 2 -2 s + 4 s yields the following zeros

ρ k = 1 ln 2 ln 1 ± i √ 7 2 + 2πik , k ∈ Z. ( 16
)
Remark 2.5. In the ∆ = 1/3 case, since all these zeros have the real part equal to 1 2 , the energy for ∆ = 1/3 exhibits only critical zeros. The last factorization we are considering in our paper corresponds to ∆ = 1/5:

E(s, 1/5) = 1 2 s+1 (3 -2 s -3 s + 6 s )ζ(s). ( 17
)
The function 3 -2 s -3 s + 6 s exhibits only off-critical zeros which can be found only numerically, e.g., s ≈ 0.635084 ± 1.07885i.

Numerical results

The starting point of our numerical determination of zeros of the energy E(s, ∆) was the case ∆ = 1/2, with the factorization form [START_REF] Hadamard | Étude sur les propriétés des fonction entières et un particulier d'une fonction considéré par Riemann[END_REF], whose spectrum of zeros involves both the critical zeros of the Riemann zeta function as well as an infinite set of off-critical zeros [START_REF] Hardy | The zeros of Riemann's zeta-function on the critical line[END_REF]. It was checked that the accuracy of determination of complex zeros by using the symbolic language Mathematica is 34-35 decimal digits for both real and imaginary components.

Then we proceeded to the left and right from this point ∆ = 1/2 by changing successively ∆ by a small amount to avoid an uncontrolled skip between neighbouring branches of zeros. Our numerical experience indicates that changing ∆ by 0.01 is certainly safe from this point of view. We observe the following:

• The zeros of E(s, ∆) form continuous non-crossing curves in the complete (ρ x , ρ y , ∆) space.

• Nevertheless, the curves may intersect in the reduced spaces (ρ y , ∆) and (ρ x , ∆), see Figures 1 and2, respectively.

Notice that the values of ∆ ∈ {1/5, 1/3, 1/2} serve as test points of our numerical calculations. The dependence of the imaginary component of zeros ρ y , in the range of its values [0, 25], on the parameter ∆ ∈ (0, 1] is pictured in Figure 1. The special cases ∆ ∈ {1/5, 1/3, 1/2}, when the energy factorizes itself onto a product of the Riemann zeta function and a simple function, are visualized by vertical dashed lines. These cases yield us precise values of the corresponding zeros and also make us sure not to miss any curve of zeros. The critical zeros with ρ x = 1 2 are denoted by red colour, all other zeros are off-critical; the off-critical zeros for the special values of ∆ ∈ {1/5, 1/2} are denoted by blue colour. We observe two types of zero curves:

• The first three "standard" zero curves ρ y (∆), which end up at the Riemann critical zeros at ∆ = 1, are represented by full symbols.

• The first three "non-standard" curves ρ y (∆) are represented by open symbols. Since in the limit ∆ → 1 -these curves tend to off-critical zeros with the divergent real component ρ x → -∞, the curves end up with crosses indicating the absence of off-critical zeros at the Riemann's ∆ = 1.

The dependence of the real part ρ x of the first three "non-standard" energy zeros on ∆ ∈ (0, 1] is presented in Figure 2 by open symbols, in close analogy with Figure 1. As before, the critical zeros with ρ x = 1 2 are denoted by red colour. We observe the following:

• For ∆ 0.75, the three curves coincide on the considered scale and go to -∞ as ∆ → 1 -.

4. Analytic results in the limit ∆ → 1 - Approaching ∆ → 1 -, with regard to Eq. ( 11) one anticipates the presence of critical zeros of the Riemann zeta function for E(s, ∆). Surprisingly, as was already indicated, there are also additional curves of off-critical zeros.

To derive coordinates of these off-critical zeros, we set ∆ = 1 -ε in (3) and expand the energy in Taylor series in the small positive ε → 0 + up to the order ε 5 .

Proposition 4.1. Let s ∈ C\{1}, then, as ε → 0 + , E(s, 1 -ε) = ζ(s) + 2 2+s -1 2 5+s s(1 + s)ζ(2 + s) ε 2 + ε 3 + 3 4 ε 4 + 1 2 ε 5 + 1 3 2 4+s -1 2 11+s s(1 + s)(2 + s)(3 + s)ζ(4 + s) ε 4 + 2ε 5 + O(ε 6 ).
Proof. It directly follows from the Taylor expansion of the Hurwitz zeta function (see e.g. [START_REF] Vepštas | An efficient algorithm for accelerating the convergence of oscillatory series, useful for computing the polylogarithm and Hurwitz zeta functions[END_REF]): for |a| < 1,

ζ(s, a) = 1 a s + ∞ n=0 (-a) n s + n -1 n ζ(s + n), (18)
where the binomial coefficient for a complex s has to be understood as

s + n -1 n = s(s + 1)(s + 2) • • • (s + n -1) n! , (19)
as well as the analytic continuation of s → E(s, ∆).

It is clear that in the limit ε → 0 + the zeros of E(s, 1-ε) coincide trivially with the critical ones of the Riemann zeta function ζ(s). Let us compute the ε → 0 + asymptotics of the other nontrivial zeros.

Theorem 4.2. The nontrivial off-critical zeros of E(s, 1 -ε) are given, as

ε → 0 + , by {ρ(k) = ρ x (k) + iρ y (k)}
+ o ε 2 ln 3 ln 2 . ( (21)
) 22
In particular, 1. Vanishing of off-critical zeros: we have lim ε→0 + ρ x (k) = -∞; 2. Asymptotic crystallization of their imaginary parts: at first order, as ε → 0 + , the imaginary parts of off-critical zeros are equidistributed on the lattice (2Z + 1)π ln 2 .

Proof. The other nontrivial zeros {ρ}, besides the one of the Riemann zeta function, correspond to solutions of the equation

2 ρ = 1 -2 2+ρ 2 5 ρ(1 + ρ)ζ(ρ + 2) ζ(ρ) ε 2 + ε 3 + 3 4 ε 4 + 1 2 ε 5 + 1 3 1 -2 4+ρ 2 11 ρ(1 + ρ)(2 + ρ)(3 + ρ)(4 + ρ)ζ(4 + ρ) ζ(ρ) ε 4 + 2ε 5 +O(ε 6 ). ( 23
)
As will be showed later, the component ρ x of ρ = ρ x + iρ y goes to -∞ as ε → 0 + . To simplify our computations, one applies the well known duality transformation

π -s/2 Γ s 2 ζ(s) = π (s-1)/2 Γ 1 -s 2 ζ(1 -s) (24)
to each Riemann zeta function in [START_REF] Spira | Zeros of Hurwitz zeta functions[END_REF]. Using then the formula Γ(x + 1) = xΓ(x), one ends up with the result

2 ρ = - 1 -2 2+ρ 2 3 π 2 ζ(-1 -ρ) ζ(1 -ρ) ε 2 + ε 3 + 3 4 ε 4 + 1 2 ε 5 + 1 3 1 -2 4+ρ 2 7 π 4 ζ(-3 -ρ) ζ(1 -ρ) ε 4 + 2ε 5 + O(ε 6 ). (25)
In the limit ε → 0 + , the r.h.s. of this equation vanishes and, consequently, the component ρ x of ρ = ρ x + iρ y must go to -∞ as indicated before. In the limit ρ x → -∞, the ratios of Riemann zeta functions in (25) can be expanded as follows

ζ(-1 -ρ) ζ(1 -ρ) = 1 + 2 ρ+1 + 3 ρ+1 + k≥4 k ρ+1 1 + 2 ρ-1 + 3 ρ-1 + k≥4 k ρ-1 = 1 + 3 2 2 ρ + 8 3 3 ρ + O(4 ρ ) (26)
and

ζ(-3 -ρ) ζ(1 -ρ) = 1 + 2 ρ+3 + 3 ρ+3 + k≥4 k ρ+3 1 + 2 ρ-1 + 3 ρ-1 + k≥4 k ρ-1 = 1 + 15 2 2 ρ + O(3 ρ ). (27)
In the leading order of the smallness parameter ε, it holds that

2 ρx+iρx = - π 2 8 ε 2 + o(ε 2 ), (28)
Since the right-hand side of this equation is real and negative, the leading order of the ρ y -component is given by 2 iρy = -1 + o(1), or, equivalently,

ρ y (k) = 1 ln 2 (2k + 1)π + o(1), k ∈ Z. (29)
This means that in the limit ε → 0 + there exists an infinite sequence of equidistant zero components along the ρ y axis. As follows from (28), the x-component of these zeros diverges logarithmically as ε → 0 + :

ρ x (k) = 2 ln 2 ln ε + -3 + 2 ln 2 ln π + O(ε). (30)
Note that the leading terms are the same for any value of k. This behavior can be seen in Figure 2. Higher orders of the expansion of ρ y (k) and ρ x (k) in ε can be obtained by inserting the leading order expressions (29) and (30) directly into the basic relation [START_REF] Titchmarsh | The Theory of The Riemann Zeta-function[END_REF]. Performing the expansion procedure in ε it is important to realize that

3 s = π 2 8 ε 2 ln 3 ln 2 exp i ln 3 ln 2 (2k + 1)π + o ε 2 ln 3 ln 2 (31)
is of order 2 ln 3/ ln 2 ≈ 3.17 > 3. After simple algebra one obtains the desired asymptotics for ρ x (k) and ρ y (k).

Comparison with our numerics. To check numerically our expansion in ε for the imaginary parts of the first three (k = 0, 1, 2) off-critical zeros, let us define the deviations from their ε = 0 values as follows

δρ y (k) := ρ y (k) - 1 ln 2 (2k + 1)π. ( 32
)
We know from ( 22) that the deviations are expected to behave in the region of the small anisotropy parameter ε → 0 + as δρ y (k) = 8 3 ln 2

Figure 1 :

 1 Figure 1: Imaginary parts ρ y of the energy zeros E(ρ, ∆) = 0. The special cases ∆ ∈ {1/5, 1/3, 1/2, 1} are visualized by vertical dashed lines. The points which correspond to critical zeros with ρ x = 1 2 are denoted by red colour.

Figure 2 :

 2 Figure 2: The dependence of the real part ρ x of zeros E(ρ, ∆) = 0, corresponding to the three "nonstandard" curves denoted by open symbols in Figure 1, on the parameter ∆ ∈ [0, 1]. The critical zeros with ρ x = 1 2 are denoted by red colour.

Figure 3 :

 3 Figure 3: Imaginary parts of the first three (k = 0, 1, 2) off-critical zeros E(ρ, ∆) = 0. The deviation δρ y (k) is defined by (32). The numerical results are depicted by open circles (k = 0), squares (k = 1) and triangles (k = 2). The plots of δρ y (k) in the region of the small anisotropy parameter ε = 1 -∆, anticipated to behave according to the asymptotic formula (33), are represented by the corresponding dashed curves.

Figure 4 :

 4 Figure 4: Real parts of the first three (k = 0, 1, 2) off-critical zeros E(ρ, ∆) = 0. The deviation δρ x (k) is defined by (34). The numerical results are depicted by open circles (k = 0), squares (k = 1) and triangles (k = 2). The plots of δρ x (k) in the region of the small anisotropy parameter ε = 1 -∆, anticipated to behave according to the asymptotic formula (35), are represented by the corresponding dashed curves.

Acknowledgements

The support received from VEGA Grant No. 2/0092/21 and Project EXSES APVV-20-0150 is acknowledged.


Laurent Bétermin

Ladislav Šamaj

Igor Travěnec

Structural transitions in interacting lattice systems

Keywords: lattice systems, double Yukawa potential, phase transition, tricritical point AMS CLASSIFICATION 74G65, 74N05, 82B26

We consider two-dimensional systems of point particles located on rectangular lattices and interacting via pairwise potentials. The goal of this paper is to investigate the phase transitions (and their nature) at fixed density for the minimal energy of such systems. The 2D rectangle lattices we consider have an elementary cell of sides a and b, the aspect ratio is defined as ∆ = b/a and the inverse particle density A = ab; therefore, the "symmetric" state with ∆ = 1 corresponds to the square lattice and the "non-symmetric" state to the rectangular lattice with ∆ = 1. For certain types of the interaction potential, by changing continuously the particle density, such lattice systems undertake at a specific value of the (inverse) particle density A * a structural transition from the symmetric to the non-symmetric state. The structural transition can be either of first order (∆ unstick from its symmetric value ∆ = 1 discontinuously) or of second order (∆ unstick from ∆ = 1 continuously); the first and second-order phase transitions are separated by the so-called tricritical point. We develop a general theory on how to determine the exact values of the transition densities and the location of the tricritical point. The general theory is applied to the double Yukawa and Yukawa-Coulomb potentials.

Introduction

Lattice systems of interacting particles are known as good models to understand physical [START_REF] Bétermin | Two-dimensional theta functions and crystallization among Bravais lattices[END_REF][START_REF] Bétermin | Lattice ground states for embedded-atom models in 2D and 3D[END_REF][START_REF] Luo | Non-hexagonal lattices from a two species interacting system[END_REF][START_REF] Luo | On minima of sum of theta functions and Mueller-Ho Conjecture[END_REF][START_REF] Sandier | From the Ginzburg-Landau model to vortex lattice problems[END_REF] or biological [START_REF] Bétermin | Theta functions and optimal lattices for a grid cells model[END_REF][START_REF] Mogilner | Mutual interactions, potentials, and individual distance in a social aggregation[END_REF] phenomena in the simplest periodic setting (see also [START_REF] Blanc | The crystallization conjecture: a review[END_REF] and references therein). In particular, the ground state approach -where a lattice energy has to be minimized in order to find the most stable state of the system -is of high interest in applied mathematics.

Considering a radially symmetric interaction potential f with parameters, a set of lattices L ⊂ R d , and defining the energy of this lattice system as

∀L ∈ L, E[L] = 1 2 p∈L\{0} f (|p|),
where | • | is the Euclidean norm on R 2 , one can ask the following important questions:

• what is the minimizer of E in the class L?

• does this minimizer changes when parameters of f vary?

• does the minimizer changes when the lattice's density varies? These questions have recently received a certain interest in the mathematical community [START_REF] Bétermin | Two-dimensional theta functions and crystallization among Bravais lattices[END_REF][START_REF] Bétermin | Optimal and non-optimal lattices for non-completely monotone interaction potentials[END_REF][START_REF] Bétermin | Lattice ground states for embedded-atom models in 2D and 3D[END_REF][START_REF] Luo | On minima of difference of Epstein zeta functions and exact solutions to Lennard-Jones lattice energy[END_REF][START_REF] Luo | On lattice hexagonal crystallization for non-monotone potentials[END_REF] as well as their more general associated problem known as crystallization (see [START_REF] Bétermin | Crystallization to the square lattice for a two-body potential[END_REF] and references therein), best packings [START_REF] Cohn | The sphere packing problem in dimension 24[END_REF][START_REF] Viazovska | The sphere packing problem in dimension 8[END_REF] and universal optimality [START_REF] Cohn | Universally optimal distribution of points on spheres[END_REF][START_REF] Cohn | Universal optimality of the E 8 and Leech lattices and interpolation formulas[END_REF]. The questions of phase transitions, i.e. studying the energy minimizers as density or parameters vary, have been shown as both interesting [START_REF] Bétermin | Local variational study of 2d lattice energies and application to Lennard-Jones type interactions[END_REF][START_REF] Travěnec | Two-dimensional Wigner crystals of classical Lennard-Jones particles[END_REF] and difficult [START_REF] Bétermin | Three-dimensional lattice ground states for Riesz and Lennard-Jones type energies[END_REF][START_REF] Luo | On minima of difference of Epstein zeta functions and exact solutions to Lennard-Jones lattice energy[END_REF] to solve.

Very few results exist in that direction. The Lennard-Jones system -i.e. for a difference of inverse power laws -is now entirely understood in dimension 2 [START_REF] Bétermin | Optimality of the triangular lattice for Lennard-Jones type lattice energies: a computer assisted method[END_REF][START_REF] Luo | On minima of difference of Epstein zeta functions and exact solutions to Lennard-Jones lattice energy[END_REF] in terms of minimizers at fixed density whereas many problems in the field are still open, with a lack of general theory. Furthermore, only few -but very important -results are known in higher dimensions, especially in dimensions 8 and 24 where best packing and universal optimality are proven [START_REF] Bétermin | Effect of periodic arrays of defects on lattice energy minimizers[END_REF][START_REF] Cohn | Universal optimality of the E 8 and Leech lattices and interpolation formulas[END_REF]. The rest of our knowledge is restricted to numerical investigations (see e.g. [START_REF] Bétermin | Three-dimensional lattice ground states for Riesz and Lennard-Jones type energies[END_REF]).

The goal of this paper is to derive technical, theoretical and numerical results for phase transition problems in lattice systems. More precisely, we are tackling here both first and second order phase transitions for smooth potentials f among two-dimensional rectangular lattices, splitting ground states into square and non-square structures. Following Landau's free energy approach in statistical physics [START_REF] Bausch | Ginzburg criterion for tricritical points[END_REF][START_REF] Landau | On the theory of phase transitions[END_REF][START_REF] Landau | Course of Theoretical Physics[END_REF][START_REF] Tolédano | The Landau theory of phase transitions[END_REF], we expand our energy E in terms of a (small) lattice parameter and we study the corresponding transition points (second-order phase transition), i.e. the value of the inverse density where transition occurs. Furthermore, we focus on important transition points called tricritical points where (continuous) second-order transition becomes of (discontinuous) first-order.

Our findings are both theoretical and numerical. After defining the classes of interaction potentials and lattices we are studying in this work, we obtain the following results:

(1) Theoretically, we show how a Taylor expansion of the energy E can be performed in terms of the lattice parameter ε (ε = 0 corresponding to the square lattice). Results are derived ensuring the existence of transition/tricritical points and showing the universal behavior of the energy minimizer in the neighborhood of these points.

(2) Numerically, we are investigating both double Yukawa and Yukawa-Coulomb potentials. Since they are highly inhomogeneous, a numerical study is performed in order to plot phase diagrams, second-order phase transition curves and tricritical points.

In both cases, asymptotics results in the neighborhood of transition/tricritical points as well as singular values of potential's parameters are derived.

The method can be generalized easily to other types of lattices, like the 2D (equilateral) triangle lattice which is the centered rectangle lattice with the aspect ratio √ 3 against the general centered rectangle lattice. The generalization of the method to 3D lattice structure is also straightforward.

The application of the theory to other types of interaction potentials, especially to Lennard-Jones potential, is of our future interest.

Plan of the paper. In Section 2, we give the precise definitions of potentials, lattices and energies we are considering. General results on phase transitions are proved in Section 3 whereas applications to double Yukawa and Yukawa-Coulomb potentials are presented in Sections 4 and 5, respectively, with both numerical and theoretical aspects.

Preliminary formalism

In this section, we briefly present the type of potential, lattices and energy we are considering in this paper.

Let us start with potentials. Our goal is to cover the main interaction potentials presented for instance in [START_REF] Kaplan | Intermolecular Interactions: Physical Picture, Computational Methods, Model Potentials[END_REF] (see also [START_REF] Bétermin | Two-dimensional theta functions and crystallization among Bravais lattices[END_REF]).

Definition 2.1 (Admissible potential). We say that f ∈ F if f : (0, ∞) → R, |f (r 2 )| = O(r -2-ε
) for some ε > 0 as r → ∞, and there exists a Radon measure µ f on (0, +∞) such that

f (r) = ∞ 0 e -r 2 t dµ f (t).
(2.1)

Furthermore, we say that f ∈ F + if µ f is non-negative.
Remarks 2.2 (Completely monotone potentials). Let us write f (r) = F (r 2 ) where F : (0, +∞) → R, then:

(1) the measure µ f is actually the inverse Laplace transform of F ;

(2) By Bernstein-Hausdorff-Widder theorem [START_REF] Bernstein | Sur les fonctions absolument monotones[END_REF], we know that f ∈ F + if and only if F is completely monotone, i.e. the derivatives of F alternate their sign: ∀k ∈ N, ∀r > 0, (-1) k F (k) (r) ≥ 0.

Examples 2.3 (Riesz and Yukawa potentials). Let us mention two important interaction potential we will consider in this work:

(1) the Riesz potential with parameter s > 0 is given by

f (r) = 1 r s , dµ f (t) = 1 Γ(s/2) t s/2-1 dt (2.2)
with Γ being the Euler Gamma function. We notice that f ∈ F + if and only if s > 2.

Nevertheless, we will also consider lower exponents by renormalizing our lattice energy (see Remark 2.10) even though f ∈ F in that case, because of its non-integrability at infinity. (2) the Yukawa potential, belonging to F + reads, for κ > 0, as

f (r) = e -κr r , dµ f (t) = 1 √ πt exp - κ 2 4t dt. (2.
3)

The set of lattice structures we are considering is defined as follows.

Definition 2.4 (Family of rectangular lattices). Let ∆ ∈ (0, 1] and A > 0. The rectangular lattice of area A and side-lengths √ A∆ and A ∆ is defined as

L ∆,A := √ A Z 1 √ ∆ , 0 ⊕ Z 0, √ ∆ ,
and its associated quadratic form is defined by

∀(j, k) ∈ Z 2 , Q ∆,A (j, k) := j √ A 1 √ ∆ , 0 + k √ A 0, √ ∆ 2 = A j 2 ∆ + k 2 ∆ ,
where

| • | is the euclidean norm on R 2 .
Remark 2.5. Notice that this family of lattices is exactly, up to isometry, the set of all rectangular lattices with fixed density. Furthermore, the particle system is invariant under rotation by the right angle π/2 which is equivalent to the interchange ∆ and 1/∆. The special case ∆ = 1 corresponds to the "symmetric" state of the square lattice, the case ∆ = 1 corresponds to the "non-symmetric" state of the rectangle lattice.

The lattice energy we are studying in this paper is the energy per point of L A,∆ interacting via potential f . Definition 2.6 (Interaction energy). Let f ∈ F, A > 0, ∆ ∈ (0, 1], then the f -energy of L A,∆ is defined by

E(A, ∆) := 1 2 ∞ ′ j,k=-∞ f Q A,∆ (j, k) = 1 2 ∞ ′ j,k=-∞ f A j 2 ∆ + k 2 ∆ . (2.4)
Remark 2.7. The prefactor 1/2 appears because each energy term is shared by two particles and the prime at the sum means that the self-energy term j = k = 0 is omitted. Note that the dependence of the energy on the parameters of the potential f is not written explicitly, for simplicity reasons.

As already done in other papers [START_REF] Bétermin | Two-dimensional theta functions and crystallization among Bravais lattices[END_REF][START_REF] Faulhuber | Optimal Gabor frame bounds for separable lattices and estimates for Jacobi theta functions[END_REF] for rectangular or general lattices, this energy can be written in terms of Jacobi theta functions thanks to the Laplace transform expression of f . Lemma 2.8 (Integral representation of the energy, see e.g. [START_REF] Bétermin | Two-dimensional theta functions and crystallization among Bravais lattices[END_REF][START_REF] Faulhuber | Optimal Gabor frame bounds for separable lattices and estimates for Jacobi theta functions[END_REF]). Let f ∈ F, A > 0, ∆ ∈ (0, 1], then we have

E(A, ∆) = 1 2 ∞ 0 θ 3 e -t/∆ θ 3 e -t∆ -1 dµ f t A , (2.5)
where

θ 3 (q) = ∞ j=-∞ q j 2 (2.6)
denotes Jacobi elliptic function with zero argument.

Remark 2.9. We use Gradshteyn-Ryzhik [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF] notation for the third Jacobi theta function. Furthermore the subtraction of -1 in the square bracket is due to the absence of the self-energy term (j, k) = (0, 0) in the sum (2.4). If µ f is absolutely continuous with respect to the Lebesgue measure, i.e. dµ f (t) = ρ f (t)dt, then we have

dµ f t A = ρ f t A dt A . (2.7)
Moreover, note that the invariance of E(A, ∆) with respect to the transform ∆ → 1/∆ is obviously ensured by the formula (2.5).

Remark 2.10. If the potential f decays to zero at large distances too slowly, the integral in (2.5) may diverge which requires a regularization (analytic continuation). Like for instance, for the Riesz potential the integral converges if s > 2. The Coulomb case with s = 1 has to be regularized by adding a uniform neutralizing background which corresponds to the subtraction of the singular term -π/t in the square bracket in the integral (2.5), as it is explained for instance in [START_REF] Travěnec | Generation of off-critical zeros for hypercubic Epstein zeta functions[END_REF].

General theory of structural transitions

To account better for the symmetry ∆ → 1/∆ of the energy (2.5), we introduce the parameter ε via ∆ = e ε .

(3.1)

Thus the ∆ → 1/∆ symmetry of the energy (2.5) is converted to the ε → -ε one of the energy

E(A, e ε ) = 1 2 ∞ 0 θ 3 e -t exp(-ε) θ 3 e -t exp(ε) -1 dµ f t A . (3.2)
At the same time, the symmetric (square lattice) value of the parameter ∆ = 1 is consistent with ε = 0.

Let us recall the universal optimality result (among lattices) concerning the square lattice, due to Montgomery [START_REF] Montgomery | Minimal theta functions[END_REF].

Proposition 3.1 (Universal optimality among rectangular lattices). If f ∈ F + , then ε = 0 is the unique minimizer of ε → E(A, e ε ) for all A > 0.
Proof. Since µ f is nonnegative and ∆ = 1 is the unique minimizer of the positive function ∆ → θ 3 (e -t∆ )θ 3 (e -t∆ -1 ) -1 as shown in [START_REF] Montgomery | Minimal theta functions[END_REF], it follows that for all ∆ > 0, E(A, ∆) ≥ E(A, 1) with equality if and only if ∆ = 1. Applying the change of variable ∆ = exp(ε) completes the proof.

The following result gives the expansion of E(A, e ε ), for fixed A and as ε → 0. Theorem 3.2 (Taylor expansion of the energy). Let f ∈ F and A > 0, then, as ε → 0,

E(A, e ε ) = E 0 (A) + E 2 (A)ε 2 + E 4 (A)ε 4 + O(ε 6 ), (3.3)
where

E 0 (A) = E(A, 1) = 1 2 ∞ 0 (θ 3 ) 2 -1 dµ f t A (3.4)
is the energy of the square lattice,

E 2 (A) = d 2 dε 2 [E(A, e ε )] |ε=0 = 1 2 ∞ 0 tθ 3 θ (1) 3 -t 2 θ (1) 3 2 + t 2 θ 3 θ (2) 3 dµ f t A , (3.5
)

and E 4 (A) = d 4 dε 4 [E(A, e ε )] |ε=0 = 1 24 ∞ 0 tθ 3 θ (1) 3 -t 2 θ (1) 3 2 + 7t 2 θ 3 θ (2) 3 + 6t 3 θ 3 θ (3) 3 -6t 3 θ (1)
3 θ

(2) 3

+t 4 θ 3 θ (4) 3 -4t 4 θ (1) 3 θ
(3)

3 + 3t 4 θ (2) 3 2 dµ f t A , (3.6)
and where the theta function and its derivatives are written for simplicity as

θ 3 := θ 3 e -t , θ (n) 3 :
= d n dt n θ 3 e -t , n ∈ N. (3.7)
Remark 3.3. Note that because of the ε → -ε symmetry of the energy, the expansion contains only even powers of ε.

Proof. We use the fact that f is analytic and f (• 2 ) is absolutely summable at infinity. Therefore, the wished expansion easily yields from the one of the theta functions product, i.e. for all t > 0 and as ε → 0, the quantity θ 3 e -t exp(-ε) θ 3 e -t exp(ε) is given by

j,k e -j 2 t exp(-ε) e -k 2 t exp(ε) = j,k e -(j 2 +k 2 )t 1 + j 2 -k 2 tε + 1 2 -j 2 -k 2 -2j 2 k 2 t + j 4 t + k 4 t tε 2 + 1 6 j 2 -k 2 1 -3j 2 t -3k 2 t -2j 2 k 2 t 2 + j 4 t 2 + k 2 t 2 tε 3 + 1 24 -j 2 -k 2 -2j 2 k 2 t + 7j 4 t + 7k 4 t -6j 6 t 2 -6k 6 t 2 +6j 4 k 2 t 2 + 6j 2 k 4 t 2 + j 8 t 3 + k 8 t 3 -4j 6 k 2 t 3 -4j 2 k 6 t 3 + 6j 4 k 4 t 3 tε 4 + O ε 6 . (3.8)
The terms of odd orders ε and ε 3 vanish because of the antisymmetry of summands with respect to the interchange of indices j and k. The terms of even orders ε 2 and ε 4 can be further simplified by using the interchange of indices j and k and the sums of type j j 2n exp(-j 2 t) with n ∈ N * are equal to (-1) n θ (n) 3 .

Remark 3.4 (Connection with Statistical Physics). The present exact expansion (3.3) has its counterpart in statistical physics where it represents a mean-field approximation for a complicated statistical model in thermal equilibrium at some temperature, known as the Landau free energy [START_REF] Bausch | Ginzburg criterion for tricritical points[END_REF][START_REF] Landau | On the theory of phase transitions[END_REF]. The parameter ε plays there the role of the order parameter which vanishes in the disordered phase and takes non-zero values in the ordered phase. The general analysis of the expansion (3.3) in the context of critical phenomena in statistical models can be found in many textbooks, see e.g. [START_REF] Landau | Course of Theoretical Physics[END_REF][START_REF] Tolédano | The Landau theory of phase transitions[END_REF]. In what follows, we shall adopt this general analysis to our ground-state problem.

By (3.3), we have that, for all A > 0, as ε → 0,

E(A, e ε ) -E(A, 1) = E 2 (A)ε 2 + E 4 (A)ε 4 + O(ε 6 ) = ε 2 E 2 (A) + E 4 (A)ε 2 + O(ε 4 ) .
Therefore, it appears that the sign of E 2 (A) determines the one of E(A, e ε )-E(A, 1) for sufficiently small values of ε > 0, giving the optimality of the square lattice when E(A, e ε ) -E(A, 1) > 0 and the optimality of a non-square one when E(A, e ε ) -E(A, 1) < 0.

Definition 3.5. Let f ∈ F\F + , then any A * such that E 2 (A * ) = 0 (3.9)
with a change of sign for E 2 at A = A * is called a transition point.

A simple condition can be written for insuring the existence of such transition point.

Proposition 3.6 (Existence of a transition point). Let f ∈ F. Then:

(1) if f ∈ F + , then there is no transition point;

(2) if f ∈ F\F + such that µ f < 0 on (0, r 0 ) or on (r 1 , +∞) for some r 0 , r 1 > 0, then a transition point exists.

Proof. It has been shown by Faulhuber and Steinerberger [START_REF] Faulhuber | Optimal Gabor frame bounds for separable lattices and estimates for Jacobi theta functions[END_REF] that ∀t > 0, tθ 3 θ

(1) 3 + t 2 θ 3 θ (2) 3 -t 2 θ (1) 3 2 > 0.
Therefore, if f ∈ F + , then µ f is nonnegative and it follows that E 2 (A) > 0 for all A > 0. Furthermore, the second point follows from [START_REF] Bétermin | Optimal and non-optimal lattices for non-completely monotone interaction potentials[END_REF]Prop. 4.4] applied to the square lattice among rectangular lattices in dimension d = 2.

Examples 3.7. Point (1) of the previous result can be illustrated by any Riesz potential, whereas point ( 2) is covered for instance by the double-Yukawa potential we study in the next section of the paper, given for all r > 0 by

f (r) = v 1 exp (-κ 1 r) r -v 2 exp (-κ 2 r) r , v 1 > v 2 > 0, κ 1 > κ 2 > 0.
We now give conditions such that a second order phase transition (see [START_REF] Šamaj | Introduction to the Statistical Physics of Integrable Many-body Systems[END_REF]), i.e. the transition for the minimizer of ε → E(A, e ε ) is continuous. Proposition 3.8 (Second order phase transition). Let f ∈ F and A * be a transition point such that:

(1) E 2 is strictly decreasing in the neighborhood of A * ;

(2) E 4 (A * ) > 0.

Then there exists A 0 > 0 and A 1 > 0 such that

• if A 0 < A < A * , then ε = 0 is the unique minimizer of ε → E(A, e ε ); • if A * < A < A 1 , a minimizer ε of ε → E(A, e ε )

satisfies the following asymptotics, as

A → A * :

ε = -dE2 dA (A * ) 2E 4 (A * ) √ A -A * + o( √ A -A * ).
Remark 3.9. The same result can be written when E 2 is strictly increasing in the neighborhood of A * , with the reverse condition E 4 (A * ) < 0, which simply exchanges the regimes of optimality.

Proof. Since a transition point exists, a simple Taylor expansion of

E 2 as A → A * reads E 2 (A) = E 2 (A * ) + dE 2 dA (A * )(A -A * ) + o(A -A * ) = - dE 2 dA (A * )(A * -A) + o(A -A * ).
By assumption, we know that dE2 dA (A * ) < 0 and let us write b := -dE2 dA (A * ) > 0 for simplicity in such a way that, as

A → A * , E 2 (A) = b(A * -A) + o(A * -A).
(3.10)

The value of ε which minimizes E(A, •) satisfies the following asymptotic stationarity condition we get from (3.3):

∂E(A, e ε ) ∂ε = 2E 2 (A)ε + 4E 4 (A)ε 3 + O(ε 5 ) = 0. (3.11)
Using (3.10), this condition, for A in the neighborhood of A * , is therefore equivalent to the couple of equations

ε = 0 or ε 2 = b 2E 4 (A * ) (A -A * ) + o(A -A * ). (3.12)
We know that the first solution is dominant in a certain region A 0 < A < A * with A 0 ≥ 0. Since ε 2 must be a real positive number and E 4 (A * ) > 0, the second solution is

ε = b 2E 4 (A * ) (A -A * ) 1/2 + o((A -A * ) 1/2 ), (3.13)
dominant in the region A * < A < A 1 for a certain A 1 > 0, which completes the proof.

Remark 3.10. Notice that indeed the transition at A = A * from ε = 0 to ε given by Eq. (3.13) is continuous. Consequently, the energy and its first derivative with respect to A are continuous as well. Furthermore, the singular behaviour ε ∝ (A -A * ) β as A → A * , A > A * , of the minimizer defines the critical exponent β which in our case acquires the mean-field value 1/2. This critical exponent is universal in the sense that it does not depend on the particular form of the interaction potential f as long as the (general) assumptions are satisfied.

Let us now assume that our interaction potential f = f α depends on a parameter α ∈ R. Thus, a transition point satisfying the assumptions of second-order transition is α-dependent and one can plot the graph of α → A * (α). We therefore get a curve of second-order transitions between square and non-square rectangular lattices as minimizers of E(A, •). In analogy with statistical mechanics, we expect only one curve of second-order transition. This curve exists say for α > α t and necessarily terminates at some α = α t once A * (α = α t ) satisfies, besides E 2 (A * (α t )) = 0, also the condition E 4 (A * (α t )) = 0 which changes sign at this point. Definition 3.11. Let f α ∈ F depending on a real parameter α. We say that A t > 0 and α t are coordinates of a tricritical point if it is a transition point satisfying both conditions

E 2 (A t (α t )) = 0 and E 4 (A t (α t )) = 0, (3.14)
with a change of sign for E 4 at A t (α t ).

Beyond this point, i.e., for α < α t , transitions between the square and rectangle lattices become of first order, i.e., there is a discontinuity in ε from 0 to a finite value. This means that our expansion (3.3) no longer applies and the curve of first-order transitions can be located only numerically by plotting the energy in the whole interval ∆ ≥ 1. Using the same method as before, we can derive the singular behavior of ε at the tricritical value of α = α t , in the region A > A t . Proposition 3.12 (Transition at the tricritical point). Let f α ∈ F depending on a real parameter α and (α t , A t ) be the coordinates of a tricritical point such that

(1) E 2 is strictly decreasing function of A in the neighborhood of A t ; (2) E 6 (A t ) = d 6 dε 6 [E(A, e ε )] |ε=0 > 0.
Therefore, there exists A 2 > 0 such that for A t < A < A 2 , a minimizer of ε → E(A, e ε ) satisfies the following asymptotics as A → A t :

ε = 4 -dE2 dA (A t ) 3E 6 (A t ) (A -A t ) 1/4 + o (A -A t ) 1/4 .
Remark 3.13. As Proposition 3.8, the same result for the A < A t regime can be written when E 2 is strictly increasing in the neighborhood of A t and E 6 (A t ) < 0.

Proof. This follows from (3.10) and (3.11) expanded to order 5 -coming from the Taylor expansion of the energy to order 6 as

E(A, e ε ) = E 0 (A) + E 2 (A)ε 2 + E 4 (A)ε 4 + E 6 (A)ε 6 + o(ε 8
)and following exactly the same lines as in the proof of Proposition 3.8.

Remark 3.14. For α = α t , we therefore get the singular behavior ε ∝ (A -A t ) 1/4 as A → A t , A > A t , so that the critical exponent at the tricritical point is β t = 1/4 which is again universal up to the assumptions.

The existence of transition and tricritical points and the second-order transition curves depends on the form of the interaction potential f (see e.g. Proposition 3.6 for the transition point). In the next two sections, we shall present explicit results for the double Yukawa and Yukawa-Coulomb potentials which are combinations of two (repulsive and attractive) terms. Lattice systems with these potentials exhibit in certain regions of model's parameters continuous as well as discontinuous structural transition from the square to rectangular lattices.

Double Yukawa potential

Potential, parameters constraints and energy. The double Yukawa interaction potential f ∈ F is defined on (0, +∞) by

f (r) = v 1 exp (-κ 1 r) r -v 2 exp (-κ 2 r) r , v 1 > v 2 > 0, κ 1 > κ 2 > 0, (4.1)
where we choose the parameters in such a way that the first repulsive term dominates at small distances whereas the second is attractive for large one. Furthermore, the corresponding measure, given by the inverse Laplace transform of r → f ( √ r), reads as

dµ f (t) = v 1 exp - κ 2 1 4t -v 2 exp - κ 2 2 4t dt √ πt . (4.2)
The potential (4.1) is supposed to have an attractive minimum f (r min ) < 0 at some r = r min . For simplicity, we set r min = 1 and f (1) = -1, reducing in this way the number of parameters four by two. We keep the independent parameters (v 1 , κ 1 ) since

κ 2 = κ 1 v 1 -e κ1 v 1 + e κ1 and v 2 = e κ2-κ1 (1 + κ 1 )v 1 1 + κ 2 . (4.3)
Since κ 2 > 0, the parameters v 1 and κ 1 are constrained to the subspace

v 1 > e κ1 κ 1 . (4.4)
The energy per particle for the double Yukawa potential reads, given (4.3), as

E(A, ∆; v 1 , κ 1 ) = 1 2 √ Aπ ∞ 0 v 1 e -κ 2 1 A 4t -v 2 e -κ 2 2 A 4t θ 3 e -t∆ θ 3 e -t/∆ -1 dt √ t . (4.5)
Notice that the exponential terms with non-zero κ 1 and κ 2 make the integral convergent for small t when θ 3 e -t∆ θ 3 e -t/∆ ∼ 1/t as t → 0.

Let us now show that the system always admits a transition point. Proof. For t > 0, we have, since v 1 > v 2 and κ 1 > κ 2 ,

v 1 exp - κ 2 1 4t -v 2 exp - κ 2 2 4t < 0 ⇐⇒ t < (κ 2 1 -κ 2 2 ) 4 ln v1 v2 =: r 0 .
It follows that µ f < 0 on (0, r 0 ) and therefore, by Proposition 3.6, E admits a transition point.

Numerical investigation.

Let us fix one of the two independent Yukawa parameters, say κ 1 = 2. By Lemma 4.1, for any given v 1 , E admits a transition point A * (v 1 ). The phase diagram in the (A, v 1 )-plane is pictured in Fig. 1. Eq. (3.9) was used to calculate the critical curve of second-order transitions defined by C := {(A * , v 1 (A * )) : A * transition point} between the square and rectangular phases, marked by open squares interconnected via a solid line. For a fixed v 1 , the square (rectangular) lattice minimizes the energy in the whole interval 0 ≤ A ≤ A * (A > A * ). We observe the following:

(1) A * → v 1 (A * ) is decreasing. (3) The critical curve C ends up at the tricritical point with the coordinates A t ≈ 2.7163619942262467 and v t 1 := v 1 (A t ) ≈ 6.7951845011079, denoted by the cross. For e κ1 /κ 1 < v 1 < v t 1 , the transition between the square and rectangular lattices is of the first order, see open circles in Fig. 1. For comparison, an artificial prolongation of the critical curve v 1 (A * ) into this region is indicated by the dashed line in the inset of the figure.

Singularity in the A > A * and A > A t regimes. To document the singular behaviour of ε ∼ ∆ -1 close to the critical curve, let us choose κ 1 = 2 and v 1 = 9.8 and the corresponding critical value of the inverse density A * ≈ 2.61449322978. For A slightly larger than A * , the rectangle energy is minimized with respect to ∆. The log-log plot of the obtained ∆ -1 versus A -A * is presented in Fig. 2. The data were fitted according to ∆ -1 ∝ (A -A * ) β . The obtained critical exponent β ≈ 0.50003 is very close to the anticipated mean-field exponent β = 1/2.

As concerns the tricritical point for the double Yukawa parameter κ 1 = 2, see the cross in Fig. 1, its coordinates A t = 2.7163619942262467 . . . and v t 1 = 6.7951845011079 . . . were calculated by using the formula (3.14). The numerical calculation of the deviation ∆ -1 in the region of the rectangular lattice A > A t close to the tricritical point was made. The log-log plot of numerical data in Fig. 3 can be fitted as ∆ -1 ∝ (A -A t ) β t , the obtained exponent β t ≈ 0.253 is reasonably close to the expected mean-field value β t = 1/4. About the transition/tricritical points and their κ 1 -dependence. The pair of coordinates of the tricritical point A t and v t 1 can be calculated by using Eq. (3.14) for any value of κ 1 ,

Yukawa-Coulomb potential

Potential, parameters constraints and energy. In the special Yukawa-Coulomb case κ 2 = 0, the potential

f (r) = v 1 exp (-κ 1 r) r -v 2 1 r (5.1)
has a minimum at r min = 1 and f (1) = -1 under conditions

v 1 = e κ1 κ 1 , v 2 = 1 + κ 1 κ 1 . (5.2)
One is left with only one independent parameter κ 1 > 0 and there are no constraints on this parameter. The associated measure is therefore

dµ f (t) = v 1 exp - κ 2 1 4t -v 2 dt √ πt . (5.3)
It is clear that f ∈ F since it is not integrable at infinity. Therefore, a regularization of the divergent lattice sum by a neutralizing background is necessary for the Coulomb term [START_REF] Travěnec | Generation of off-critical zeros for hypercubic Epstein zeta functions[END_REF]. The In Fig. 8, ∆ -1 is plotted as the function of (A -A t ) in order to check the tricritical exponent in the Yukawa-Coulomb case. The numerical data are represented by open circles. Fitting the log-log plot via ∆ -1 ∝ (A -A t ) β t , represented by the dashed line, yields the tricritical exponent β t = 0.2498 which is very close to the anticipated mean-field one β t = 1/4.

Figure 1 .

 1 Figure 1. Phase diagram in the (A, v 1 )-plane for the double Yukawa potential with the fixed parameter κ 1 = 2. The critical curve of second-order transitions between the square and rectangular lattices is marked by open squares interconnected via a solid line. The cross marks the tricritical point. The discontinuous first-order transitions are indicated by open circles. The dashed line in the inset, given by E 2 (A) = 0, represents an artificial extension of the second-order transitions to the region dominated by first-order transitions. The dash-dotted straight line v 1 = e κ 1 /κ 1 = e 2 /2 denotes the border of v 1 -values.

Lemma 4 . 1 (

 41 Existence of transition points). For any fixed (v 1 , κ 1 ), E admits a transition point A * (v 1 , κ 1 ).

Figure 2 .

 2 Figure 2. The double Yukawa potential with the parameters κ 1 = 2 and v 1 = 9.8, corresponding to the critical value of A * ≈ 2.61449322978. The log-log plot of ∆ -1 versus A -A * is presented, numerical data are represented by open circles interconnected via a dashed line. The fitting via ∆ -1 ∝ (A -A * ) β yields the critical exponent β * ≈ 0.50003.

( 2 )

 2 There exists a finite "minimal" value of A * , A * min ≈ 2.18626, such that lim A→A * min A>A * min v 1 (A * min ) = +∞.

Figure 3 .

 3 Figure 3. The double Yukawa potential with the parameter κ 1 = 2, the corresponding tricritical point has coordinates A t = 2.7163619942262467 . . . and v t 1 = 6.7951845011079 . . .. The log-log plot of ∆ -1 versus A -A t is presented in the rectangular region A > A t , numerical data are denoted by open circles interconnected via a dashed line. The fitting via ∆ -1 ∝ (A -A * ) β yields the tricritical exponent β t ≈ 0.253.

Figure 4 .

 4 Figure 4. The components (A t , v t 1 ) of the tricritical points as the functions of the Yukawa parameter κ 1 . The tricritical point exists only if κ 1∈ (κ L 1 , κ U 1 )where κ L 1 ≈ 1.436 and κ U 1 ≈ 2.03414 are the lower and upper limits, respectively. The values of A t are represented by squares and the ones of v t 1 by circles. The dash-and-dot line corresponds to the curve v 1 = e κ 1 /κ 1 , there are no physical solutions below it.

1 *Figure 5 .

 15 Figure 5. The dependence of the minimal value of the critical inverse density A * min , at which a second-order transition takes place in the limit of v 1 → ∞, on the Yukawa parameter κ 1 > 0. The limiting values of A * min (κ 1 ) are: A * min (κ 1 → ∞) = 1 and A * min (κ 1 → 0 + ) = 5.71344 . . ..

Figure 7 .

 7 Figure 7. Phase diagram of the Yukawa-Coulomb model close to the tricritical point denoted by the cross. The critical line of second-order transitions is denoted by open squares, first-order transitions are indicated by open circles. The dashed line corresponds to an artificial prolongation of the critical line of second-order transitions to the region where first-order transitions are dominant.

1 ∆Figure 8 .

 18 Figure 8. The plot of the deviation ∆ -1 versus A -A t for the Yukawa-Coulomb model close to the tricritical point.

Acknowledgment

The support received from the project EXSES APVV-20-0150 and VEGA Grant No. 2/0092/21 and is acknowledged.

the obtained results are presented in Fig. 4. It turns out that a physically meaningful solution for the tricritical point exists only if κ 1 ∈ (κ L 1 , κ U 1 ) where κ L 1 ≈ 1.436 and κ U 1 ≈ 2.03414 are the lower and upper limits, respectively. While v t 1 is a decreasing function of κ 1 , A t is a slowly increasing function of κ 1 . The value of v t 1 diverges when κ 1 approaches its lower limit κ L 1 from above. Consequently, all transitions are of first order (discontinuous) for 0 < κ 1 ≤ κ L 1 . The upper limit κ U 1 is given by the intersection of the v t 1 -curve with the dash-and-dot line v 1 = e κ1 /κ 1 which is the border of the accessible subspace for Yukawa parameters (4.4). All transitions are of second order (continuous) for κ 1 > κ U 1 ; for every κ 1 from this interval, the critical line goes down up to the border line v 1 = e κ1 /κ 1 , leaving no space for a tricritical point and the corresponding first-order transitions.

Recall that for κ 1 = 2 the "minimal" critical value of A * at which v 1 (A * ) goes to infinity was found numerically to be nonzero, in particular A * min ≈ 2.18626. In what follows, we aim at deriving an exact formula for A * min (κ 1 ) for any value of κ 1 > 0. According to (4.3), for a fixed value of κ 1 and in the limit of large v 1 , the Yukawa parameter κ 2 behaves as

and the Yukawa parameter v 2 is given by

Since

using the asymptotic expansions (4.6) and (4.7), the expression into brackets can be expressed as

With regard to Eq. (3.5), the critical condition E 2 (A * ) = 0 takes in the limit

This equation determines for any κ 1 > 0 the exact value of the critical inverse density A * min (κ 1 ) at which v 1 → ∞. In particular, at κ 1 = 2 the obtained A * min ≈ 2.186262818188 agrees with the previous numerical estimate A * min ≈ 2.18626. The monotonous decay of A * min with increasing the Yukawa parameter κ 1 > 0 is presented in Fig. 5

by the solid line connecting data (open circles).

The function A * min (κ 1 ) tends to unity in the limit κ 1 → ∞. In the opposite limit

3 -t θ 1) versus ∆ for the Yukawa parameter κ 1 = 2.0365, just below its tricritical value κ t 1 = 2.036517758847 . . . where a first-order transition between the square and rectangular lattices exists. Three values of the inverse particle density A are considered: A = 2.79544356250 when the square lattice with ∆ = 1 prevails, the first-order transition value Atrans = 2.795443562576 where the energies of the square and rectangular lattices, separated by energy barriers, coincide and A = 2.795443562606 where the rectangular case is dominant via a jump in ∆.

energy then reads as

The neutralizing background manifests itself as the addition of the singular term -π/t in the last square bracket, which ensures the convergence of the integral. It is again straightforward to show that transitions point exists by direct application of Proposition 3.6 in the same way that we did in Lemma 4.1 for the κ 2 > 0 case.

Numerical investigation. Solving the couple of equations (3.14) with κ 2 = 0, v 1 and v 2 given by (5.2), one gets the tricritical point for the Yukawa-Coulomb interaction: A t ≈ 2.795433950879 and κ t 1 ≈ 2.036517758847. Let us make a small step into the region κ 1 < κ t 1 , where a step-wise first-order transition exists, say κ 1 = 2.0365. We observe the following:

(1) As is seen in Fig. 6, there is just one energy minimum at ∆ = 1 for A slightly below the first-order transition value A trans = 2.795443562576 and the square lattice prevails. (2) For A slightly larger than A trans , one gets two equivalent minima with ∆ > 1, related via the symmetry ∆ → 1/∆, meaning that the rectangular case is dominant via a jump in ∆. (3) Exactly at A trans we have three equivalent minima, i.e., the energies of the square and rectangular lattices, separated by energy barriers, coincide.

The phase diagram of the Yukawa-Coulomb model is pictured in Fig. 7. The second-order


Bernard Diény
email: bernard.dieny@cea.fr

Olivier Fruchart

Ernesto E Marinero

B Dieny

Comparison of stable spin textures in in-plane vs. out-of-plane magnetized exchange-biased multilayers

Introduction:

The shift of the hysteresis loop of a ferromagnet (FM) when grown adjacent to an antiferromagnet (AFM), the exchange bias effect, was discovered by Meiklejohn and Bean in 1957 [START_REF] Meiklejohn | New magnetic anisotropy[END_REF].

It arises from interfacial magnetic exchange coupling between the constituent magnetic layers. These thin-film FM/AFM bilayers are model systems to investigate and understand the underlying physics of exchange bias. The unidirectional pinning of the ferromagnet resulting from exchange bias can be established in a FM/AFM bilayer by cooling it in an external magnetic field from above the AFM blocking temperature or by applying a magnetic field during thin film growth [START_REF] Meiklejohn | New magnetic anisotropy[END_REF]. As a result, the hysteresis loop of the FM layer exhibits a field-shift away from the origin, whose magnitude is called the exchange bias field. This effect has been of scientific interest for decades, and it depends on the magnetic properties and the microstructure of the FM and AFM materials. It is also exploited in a wide variety of magnetic devices to pin the magnetization of the reference layer in metallic spin-valves or in-plane magnetized magnetic tunnel junctions (MTJ) used in particular as hard-disk-drive (HDD) readers [START_REF] Fullerton | Spintronics, Magnetoresistive Heads, and the Emergence of the Digital World[END_REF], as well as magnetic field sensors [START_REF] Freitas | spintronic sensors[END_REF], amongst others. Besides, the electrical manipulation of exchange bias opens new possibilities for the development of low-energy consumption spintronic devices [START_REF] Allibe | [END_REF].

Numerous studies have investigated the origins and engineering of exchange bias in FM/AFM structures. Salient findings from these studies include: exchange coupling between the ferromagnet and the antiferromagnet is long-ranged in nature (it has been observed across 5.5 nm-thick nonmagnetic spacers) [START_REF] Gökemeijer | Long-range exchange coupling between a ferromagnet and an antiferromagnet across a nonmagnetic spacer layer[END_REF]; during magnetization reversal, the AFM spins may form exchange-spring domain wall structures extending along the thickness on both sides across the FM/AFM interface; the external field was found to move only the ferromagnetic part of the hybrid domain walls, while the antiferromagnetic parts of the hybrid domain walls defining the AFM domains remain unaffected [START_REF] Chien | Hybrid domain walls and antiferromagnetic domains in exchange-coupled ferromagnetic/antiferromagnetic bilayers[END_REF]; 360° domain-wall rings may form upon field-driven magnetization reversal of in-plane-magnetized films, resulting from spatial fluctuations of the exchange-bias field [START_REF] Cho | Characteristics of 360°-domain walls observed by magnetic force microscope in exchange-biased NiFe films[END_REF]. These rings constitute a type of magnetic texture topologically protected by its close loop shape and remaining stable up to fields of several tens of mT.

Although exchange bias has been mainly studied and used in spin-valves and tunnel junctions magnetized in-plane, it is not limited to such systems. The existence of perpendicular exchange biasing for FM/AFM multilayers in Co/Pt multilayers exchanged coupled to CoO thin layers was demonstrated by Maat et al. [START_REF] Maat | Perpendicular exchange bias of Co/Pt multilayers[END_REF]. The interfacial exchange energies were found to be twice larger for in-plane vs. outof-plane field cooling geometry. This was attributed to differences in spin direction in the AFM after cooling. Subsequently, perpendicular exchange bias was reported with FeMn [START_REF] Garcia | Exchange bias in (Pt/CoFe)/FeMn multilayers with perpendicular magnetic anisotropy[END_REF] as well as in synthetic (Co/Pt)/Ru/(Co/Pt)/FeMn exchange biased ferrimagnetic pinned layer [START_REF] Sort | Pinned synthetic ferrimagnets with perpendicular anisotropy and tunable exchange bias[END_REF].

Exchange bias was also observed in systems consisting of layers with orthogonal anisotropies.

Marinero and coworkers studied the magnetic exchange interactions in Fe/TbFe wherein the Fe layer has in-plane magnetization (Fe) while the ferromagnetic TbFe exhibits out-of-plane anisotropy [START_REF] Marinero | Magnetic interactions in Fe/TbFe bilayers[END_REF][START_REF] Aeschlimann | Exchange coupling of contacted ferromagnetic films: Fe on amorphous TbFe[END_REF].

The interfacial coupling was shown to result in canting of the TbFe spins away from the perpendicular easy axis (perpendicular), and be responsible for the observed exchange bias effect [START_REF] Siegmann | Origins of surface induced magnetic structures[END_REF]. Orthogonal exchange bias was also reported in NiFe magnetized in-plane and exchange-coupled to out-of-plane magnetized (Pt/Co) multilayers [START_REF] Bollero | Exchange bias in annealed [Pt/Co]/NiFe systems[END_REF][START_REF] Bollero | Origin and control of bias-like phenomenon in coupled ferromagnetic [Pt/Co]/NiFe bilayers[END_REF]. In these systems, exchange bias appears when the (Pt/Co) multilayer is in a multidomain configuration and results from interfacial exchange interaction between the in-plane polarized closure domains of the (Pt/Co) multilayer and the NiFe magnetization. Similarly, exchange bias was observed in FeCoB in-plane magnetized layer, exchange coupled to out-of-plane magnetized SmCo5 ferromagnet [START_REF] Bollero | An extraordinary chiral exchange-bias phenomenon: engineering the sign of the bias field in orthogonal bilayers by a magnetically switchable response mechanism[END_REF].

In recent years, the concept of topological protection of spin textures has received great interest, often involving closed shapes of magnetization vector fields, e.g. on spheroids in 3D or along loops in 2D. This is in particular the case of skyrmions, found in systems lacking inversion symmetry and displaying Dzyaloshinskii-Moriya interactions [START_REF] Moriya | Anisotropic superexchange interaction and weak Ferromagnetism[END_REF][START_REF] Dzyaloshinskii | Thermodynamic theory of weak ferromagnetism in antiferromagnetic substances[END_REF]. These spin textures take the form of tubes in bulk systems [START_REF] Rößler | Spontaneous skyrmion ground states in magnetic metals[END_REF][START_REF] Mühlbauer | Skyrmion lattice in a chiral magnet[END_REF][START_REF] Pappas | Chiral paramagnetic skyrmion-like phase in MnSi[END_REF][START_REF] Yu | Real-space observation of a two-dimensional skyrmion crystal[END_REF] and bubbles in two-dimensional multilayered systems [START_REF] Heinze | Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions[END_REF][START_REF] Boulle | Room-temperature chiral magnetic skyrmions in ultrathin magnetic nanostructures[END_REF]. They were proposed to encode information at high density [START_REF] Fert | Skyrmions on the track[END_REF]. However, despite also displaying chiral features, the origin and magnetic stability of the 360° domain wall rings in in-plane magnetized exchange biased systems [START_REF] Fert | Skyrmions on the track[END_REF][START_REF] Schafer | Domain and Domain Wall Observations in Sputtered Exchange-Biased Wedges[END_REF][START_REF] Zhao | Study of 360° domain walls in NiFe/NiO film by tip-sample interaction on magnetic force microscope[END_REF][START_REF] Chapman | Magnetization reversal process in TbCo-biased spin valves[END_REF][START_REF] King | Magnetization reversal of NiFe films exchange-biased by IrMn and FeMn[END_REF][START_REF] Boeve | Creation and annihilation of 360 domain walls in magnetic tunnel junctions with exchange-biased artificial ferrimagnet[END_REF][START_REF] Kovács | Reversal mechanism of CoFeB/IrMn bilayers observed by Lorentz electron microscopy[END_REF][START_REF] Dean | The formation mechanism of 360 domain walls in exchanged-biased polycrystalline ferromagnetic films[END_REF] derive from different physical mechanisms and have a different topology (Fig. 1a). In this context of similarities yet with key differences, it is interesting to note that skyrmions have been stabilized at room temperature and zero applied magnetic field thanks to the use of exchange bias [START_REF] Rana | Room-temperature skyrmions at zero field in exchange-biased ultrathin films[END_REF] (Fig. 1b). 3) (image adapted from Ref [START_REF] Rana | Room-temperature skyrmions at zero field in exchange-biased ultrathin films[END_REF]).

a) b) 400nm 1mm

In this article, we discuss and compare spin textures formed in exchange biased magnetic structures with in-plane and out-of-plane magnetic anisotropy.

360 Domain Wall (360 DWs) segments and rings in exchange biased systems with in-plane magnetization:

The formation of 360 domain walls in exchange-biased FM films have been a subject of interest over the last 30 years for both their practical implications in spintronic devices as well as to understand the fundamental physical mechanisms and material properties that control their formation, shape, size, and stability against temperature, applied magnetic fields and spin currents. Indeed, while 360° DWs in the form of segments and bubbles have been observed in ferromagnetic single films [START_REF] Smith | Noncoherent Switching in Permalloy Films[END_REF][START_REF] Wade | Some factors in the easy axis magnetization of permalloy films[END_REF][START_REF] Masseboeuf | Probing magnetic singularities during magnetization process in FePd films[END_REF] and strips [START_REF] Zhang | 360° domain walls: stability, magnetic field and electric current effects[END_REF], they are much more common in exchange-biased films. In these systems, 360 DWs cause noise in spin-valve magnetic sensors as they constitute micromagnetic objects with strong irreversibility. Similarly, in MTJs, they are known to cause degradation in signal to noise ratio [START_REF] Dzyaloshinskii | Thermodynamic theory of weak ferromagnetism in antiferromagnetic substances[END_REF]. 360° DWs can be observed in the form of open loop segments (see for instance at the center of Fig. 1a) or rings. In the case of segments, the two ends are pinned by defects, while the 360° DW rings are stabilized by their close loop shape combined with the easy-plane shape anisotropy. The rings are usually more stable than the segments. Once such ring structures are formed, they can remain stable at applied magnetic fields up to several tens of mT, hence hindering the operating specification range of sensor devices [START_REF] Fert | Skyrmions on the track[END_REF].

Dean et al. studied the formation mechanism of 360 DWs in an exchange-biased bilayer of CoFeB/IrMn [START_REF] Cho | Characteristics of 360°-domain walls observed by magnetic force microscope in exchange-biased NiFe films[END_REF]. They incorporated the measured magnetic properties of the bilayer into a micromagnetic model of exchange bias for granular AFM films and results were compared to imaging of the micro-magnetic structure determined from Lorentz TEM. The following are the salient conclusions of their study: i) to generate a stable 360° DW, one requires multiple nuclei of reversed domains having rotated all clockwise or anticlockwise, generating 180 Néel domain walls of same sense of rotation facing each other in the various nuclei, which under a sweeping magnetic field are pushed together to form a 360° DW. Application of the magnetic field at a small angle with respect to the EB direction promotes the formation of nuclei with the same sense of rotation; ii) Multiple nucleation regions emerge naturally as the exchange field acting on the FM is spatially non-uniform. This is a consequence of the distribution of anisotropy axis in the IrMn grains. In contrast, samples with minimal crystalline in-plane dispersion, such as epitaxial Fe (001)[100]//IrMn(001) [110], bilayers are less prone to form nucleation regions with reversed magnetization, and experiments confirm that no 360° DW are observed in such samples [START_REF] Wang | Exchange bias in epitaxial Fe/Ir0.2Mn0.8 bilayers grown on MgO (0 0 1)[END_REF]; iii) the stability and density of 360° DW arises from the balance between intergrain exchange and anisotropy energies.

The internal spin-structure of the ring could be determined from Lorentz imaging as well as MFM imaging, as illustrated in Fig. 2 [START_REF] Gillies | Magnetization reversal mechanisms in NiFe/Cu/NiFe/FeMn spin-valve structures[END_REF]. The domain wall ring consists of a closed loop 360° Neel wall, i.e., magnetization rotates with a unique sense by 360° from the inner edge to the outer edge of the wall. However, at a given radial position in the ring, the spin orientation is uniform around the ring. This means, 360° domain walls located diametrically opposed in the ring have opposite chirality, contrary to the case of skyrmions. The magnetostatic charges generated by this magnetization configuration are represented in Fig. 2 as well as the magnetic contrast expected from this charge distribution when imaging this spin texture with an MFM tip magnetized out-of-plane, which fits the experimental result shown in Fig. 1a. The charge distribution is therefore chiral, which again is contrary to the case of skyrmions, whose charge distribution is invariant upon rotation. Such 360° domain wall ring textures in in-plane exchange-biased films are quite stable, requiring fields of several tens of mT to be annihilated. However, they are not topologically protected due to the uniform magnetization for any given radial distance. Expressed in topological terms, the winding number [START_REF] Braun | Topological effects in nanomagnetism: from superparamagnetism to chiral quantum solitons[END_REF] 𝑆 = 1 4𝜋 ∫ 𝑛(𝑥, 𝑦) d𝑥d𝑦 with 𝑛(𝑥, 𝑦) = 𝑚 • (𝜕 𝑥 𝑚 × 𝜕 𝑦 𝑚) is equal to zero since these spin textures are fully in-plane. This spin texture could in principle unwind through the 360° rotation of the central part. In practice, this does not occur as, at zero field, the core magnetization as well as the edge magnetization are maintained along the exchange bias field direction. However, this spin texture can also unwind by out-of-plane 180° flipping of the magnetization located at mid-width of the 360° ring.

The energy barrier protecting these rings therefore arises from the in-plane confinement of the spins, and not from the creation of a Bloch point micromagnetic defect, as in the case of so-called topological protection.

Fig.2: magnetization and charge structure 360° DW ring. The local density of magnetostatic charges determined by the divergence of magnetization is represented in red (positive charges) and green (negative charges). The white and grey contrasts indicate the expected MFM contrasts corresponding

to those experimentally observed in Fig. 1a.

Note that in a given experimental realization in a sample, all rings have the same left hand/right hand dark/bright contrast (See for instance in Fig. 1a, this is also evident in similar images reported in literature [START_REF] Cho | Characteristics of 360°-domain walls observed by magnetic force microscope in exchange-biased NiFe films[END_REF][START_REF] Zhao | Study of 360° domain walls in NiFe/NiO film by tip-sample interaction on magnetic force microscope[END_REF][START_REF] King | Magnetization reversal of NiFe films exchange-biased by IrMn and FeMn[END_REF][START_REF] Kovács | Reversal mechanism of CoFeB/IrMn bilayers observed by Lorentz electron microscopy[END_REF][START_REF] Dean | The formation mechanism of 360 domain walls in exchanged-biased polycrystalline ferromagnetic films[END_REF]). This implies that all 360° DW have the same sense of rotation from inside to outside the rings. This means that all (or at least a majority of) nuclei must have rotated along the same sense at the nucleation stage; otherwise there would be an equal share of right and left contrast spirals. Note that, contrary to the case of skyrmions in perpendicularly-magnetized films, the Dzyaloshinskii-Moriya interaction cannot explain this chirality in the charge distribution. Indeed, as the DMI vectors are oriented in the plane of the films for interfacial-based DMI, there is no DMI energy for purely in-plane-magnetized films.

As already mentioned, the selection of a given chirality of the rings is generally understood by considering their formation mechanism, as there always exists a slight misalignment between the exchange anisotropy field and the applied field, typically of a few degrees. The exchange anisotropy field is commonly set by the magnetic field applied during cooling the AFM/FM system from above the AFM blocking temperature [START_REF] Meiklejohn | New magnetic anisotropy[END_REF]. Let's then assume that an hysteresis loop measurement is performed with the field applied slightly tilted by a few degrees away from the field cooling direction in the anticlockwise direction. Starting from positive saturation and decreasing the field, the magnetization will tend to rotate towards the exchange anisotropy field, i.e. in the clockwise direction. In the ideal situation of a perfectly uniform exchange field, all magnetization vectors would rotate along the same direction, leading to a coherent reversal or nucleation-propagation reversal. However, it is known that the exchange anisotropy field is inhomogeneous due to the distribution of the anisotropy axis in the AFM grains, on account for example of variations in crystalline growth orientation, and to the frustration of exchange interactions at the FM/AFM interface. This was emphasized by Malozemoff, who proposed a random field model of exchange anisotropy at a rough FM/AFM interface [START_REF] Malozemoff | Random-field model of exchange anisotropy at rough ferromagneticantiferromagnetic interfaces[END_REF]. This random field adds up to the average exchange anisotropy field. Wherever the sum of this random field plus the exchange anisotropy field makes an opposite angle with respect to the applied field direction, a correlated area of magnetization forms whose moment rotates in the direction opposite to that of the majority of other moments, say, counterclockwise. As the field is progressively reversed, both the net magnetization and the correlated areas rotate away from their initial direction. Due to the misalignment of the applied field with the mean exchange anisotropy field, one of the two subsets of correlated areas forms a majority, while the other is a minority, thus progressively forming islands inside a matrix deriving from the former correlated areas. When the field gets close to negative saturation, each nucleated island is surrounded by a closed loop 360° Neel DW. These DWs have all the same left hand/right hand chirality, which is determined by the sign of the angular misalignment between the net exchange bias field direction and the applied field. The fact that most of these islands have the shape of circular rings of approximately the same diameter (see Fig. 1a) may originate from an energy minimization, the shrinking being stopped at some point due the inability to unwind associated with the in-plane spin confinement associated with the easy-in plane shape anisotropy, or curvature-induced magnetostatic repulsion.

Skyrmions and skyrmoniums in exchange biased FM with perpendicular anisotropy:

Skyrmions have the same topology as a bubble with 180° wall in a perpendicularly-magnetized medium, with rotation symmetry. In thin films, the domain wall is of Néel type despite their dipolar energy, promoted by the Dzyaloshinskii-Moriya interaction (DMI). DMI makes skyrmions a chiral object, in which Néel walls point radially outwards or inwards depending on the sign of the DMI interaction, and whether the skyrmion has a center pointing up surrounded by down, or the reverse [START_REF] Thiaville | Dynamics of Dzyaloshinskii domain walls in ultrathin magnetic films[END_REF][START_REF] Chen | Novel chiral magnetic domain wall structure in Fe/Ni/Cu(001) films[END_REF]. Skyrmions are characterized the winding number +1, reflecting their chiral nature and the fact that magnetization spans 4𝜋 on the unit sphere over the entire spin texture.

Isolated skyrmions are metastable objects, which may be obtained following specific magnetization or current processes. A large DMI interaction is desirable to promote the stability of skyrmions, especially under dynamic conditions. Above a given threshold of DMI interactions, the ground state is stripe domains under zero magnetic field, a periodic network of skyrmions for moderate field, and isolated metastable skyrmions shortly before saturation. Skyrmions are therefore quite sensitive to an applied magnetic field, which influences their diameter through compression, similar to the 360° DW rings in in-plane EB systems. Exchange bias of FM with perpendicular anisotropy was recently utilized [START_REF] Rana | Room-temperature skyrmions at zero field in exchange-biased ultrathin films[END_REF][START_REF] Chen | Room-temperature skyrmion ground state stabilized through interlayer exchange coupling[END_REF] to stabilize skyrmions at room temperature and zero applied external magnetic fields (Fig. 1b). The interlayer exchange coupling acts as an effective magnetic field that breaks the symmetry of "up" and "down" domains, resulting in the formation of skyrmions inside one single domain with opposite magnetization. In this study, IrMn was used as the AF, which is a very common choice. Ma et al.

calculated and observed that Dzyaloshinskii-Moriya Interactions (DMI) exist across FM/AFM interfaces of IrMn/CoFeB/MgO multilayer thin films [START_REF] Ma | Dzyaloshinskii-Moriya Interaction across an Antiferromagnet-Ferromagnet Interface[END_REF], with a DMI vector in the plane of the interface. This extra DMI adds up to the usual one of skyrmionic films promoted by other interfaces, and contributes to impose a uniform chirality of the domain walls as in standard skyrmions. In particular, note that the contrast is rotationally-invariant, contrary to 360° domain-wall rings for in-plane magnetized exchange-biased systems (fig. 2). Skyrmionic textures with 360° domain walls also exist, called skyrmoniums (see fig. 3). These are again metastable ring textures, now spanning 8𝜋 over the unit sphere for magnetization and winding number 2, which can be promoted by specific magnetizing or current history [START_REF] Zhang | Real-Space Observation of Skyrmionium in a Ferromagnet-Magnetic Topological Insulator Heterostructure[END_REF]. Their distribution of dipolar charges is rotational-invariant, similar to skyrmions, and therefore distinct from the in-plane 360° DW rings in EB systems.

Conclusion:

We have compared spin textures that develop in exchange biased systems either in-plane or outof-plane magnetic anisotropy. For in-plane exchange biased systems, 360° Neel domain wall rings may form, which are quite stable when subjected to an external magnetic field. Their stability originates from the closed loop shape of the domain wall combined with the in-plane confinement of the spins due to the easy-plane shape anisotropy. Their spin texture is not chiral and has zero winding number, but the distribution of associated dipole charges is chiral, reflecting the direction of rotation of the inplane magnetization during the initial stages of the reversal. Out-of-plane exchanged-biased systems may be influenced by Dzyaloshinskii-Moriya Interactions, whereas DMI are absent in systems with inplane anisotropy. As a result, DMI promotes the formation of skyrmions and skyrmoniums, with 180° and 360° DWs and winding number 1 and 2, respectively, similar to thin films with no EB. In contrast to in-plane magnetized systems, their distribution of dipolar charges is not chiral.

Fig. 1 :

 1 Fig.1: MFM images of (a) 360° Néel-domain-wall rings in Co(3)/IrMn(7), displaying a chiral charge distribution; (b) skyrmions stabilized at room temperature and at zero field in Ta/Pt(3)/Co(0.3)/Ni80Fe20(2.1)/Ir20Mn80(4.2 )/Pt(3) (image adapted from Ref[START_REF] Rana | Room-temperature skyrmions at zero field in exchange-biased ultrathin films[END_REF]).

Fig. 3 :

 3 Fig.3: Skyrmonium structure. From[START_REF] Zhang | Real-Space Observation of Skyrmionium in a Ferromagnet-Magnetic Topological Insulator Heterostructure[END_REF]

Acknowledgements: E.M. acknowledges financial support from CNRS for a sabbatical at SPINTEC. We acknowledge fruitful discussions with our colleague O. Boulle.

Data availability:

The authors make the data supporting this paper available upon reasonable request.

Conflict of interest:

The authors that they have no conflict of interest in relation to this research.