Keywords: Hospital, random arrivals, random exits, types of patients, available care time, nurses working conditions, Poisson distribution. Subject Classification: 60G50, 91A40, 91B32, 91B70

Introduction

As part of his work on consumption, poverty, and well-being, which culminated in the award of the Sveriges Riksbank Prize in Economic Sciences in Memory of Alfred Nobel in 2015, Angus Deaton wrote a book entitled "The Great Escape", providing an overview of developments in humanity over the last two hundred and fty years. Referring explicitly to the title of John Sturges' 1963 lm, Deaton describes how humanity has escaped hardship and premature death over this period: "Life is better now than at almost any time in history. More people are richer and fewer people live in dire poverty. Lives are longer and parents no longer routinely watch a quarter of their children die. Yet millions still experience the horrors of destitution and of premature death. The world is hugely unequal." [2013, p. 9]. This spectacular evolution, like the glaring inequalities we see between countries and populations, naturally lies in the evolution of various components but, as Deaton points out, income and health are two of the main components of human well-being. Moreover, in the labor-intensive healthcare sector, there is a correlation between the number of doctors and patient prognosis, as demonstrated, for example, in the study by [START_REF] Needleman | Nurse-stang levels and the quality of care in hospitals[END_REF]. The role of nurses is particularly crucial since, as stated by the World Health Organization (WHO [2022]), through their work they promote health, prevent disease, provide primary and community care, deliver care in emergency situations, and their participation is essential to achieving universal health coverage. Nurses are thus at the heart of any healthcare system.

Unfortunately, there is a worldwide shortage of healthcare personnel, and the COVID-19 crisis has exacerbated these shortages. This is particularly true for nurses and midwives, who account for over 50% of the world's current unmet need for healthcare personnel (there is currently a global shortage of 900,000 midwives and 6 million nurses); although the situation is highly uneven, with almost 90% of these shortages concentrated in low-income countries, particularly in Africa and South-East Asia, as the United Nations Regional Information Center (UNRIC [2021]) notes. Looking ahead, WHO estimates that an additional 9 million nurses and midwives will need to be recruited worldwide by 2030 (WHO [2022]) and, [START_REF] Scheer | Projecting shortages and surpluses of doctors and nurses in OECD[END_REF] have predicted a shortfall of 2.5 million nurses in 2030 for 23 OECD countries. In short, demand exceeds supply everywhere.

Those shortages have two major implications.

First, they call into question decades of improvement in the provision of care and patient management, by increasing the risk of death and increasing the loss of opportunity (e.g., [START_REF] Twigg | The impact of undestaed shifts on nurse-sensitive outomes[END_REF], [START_REF] Griths | Registered nurse, healthcare support worker, medical stang levels and mortality in English hospital trusts: a cross-sectional study[END_REF], Haegdorens et al. [2019], [START_REF] Needleman | Association of registered nurse and nursing support stang with inpatient hospital mortality[END_REF], Keck School of Medicine of USC [2023]). By way of illustration, the Royal College of Emergency Medicine (RCEM [2023]) estimates the number of deaths in England at 442 per week for the year 2022 due to deciencies in emergency care. Equally worrying, on January 8, 2023, the French Minister of Health declared in a televised interview: "We have fewer professionals, fewer nurses, fewer orderlies, so we have beds that are closed. We are not closing them for the sake of closing them, we are closing them because we have fewer sta" (France Info [2023]). And although sta shortages alone do not explain all the hospital bed closures,2 they do contribute to them.

Second, through the degradation of their working conditions, those shortages also have negative impacts on the physical, emotional, and mental health of caregivers (sleep disorders, stress, ethical suering, burnout, depression, injuries, etc.) [START_REF] Rotenstein | Prevalence of Depression, Depressive Symptoms, and Suicidal Ideation Among Medical Students: A Systematic Review and Meta-Analysis[END_REF], [START_REF] Duarte | Male and Female Physician Suicidality: A Systematic Review and Meta-analysis[END_REF], [START_REF] Hardy | Comparison of burnout, anxiety and depressive syndromes in hospital psychiatrists and other physicians: Results from the ESTEM study[END_REF], [START_REF] Sexton | Emotional Exhaustion Among US Health Care Workers Before and During the COVID-19 Pandemic, 2019-2021[END_REF]). These elements thus contribute to increased absenteeism, sick leave, and even abandonment of their profession by many nurses, as highlighted, for example in the USA by the National Council of State Boards of Nursing (NCSBN [2023]), which reports that 100,000 graduate nurses have left the profession since 2020, and that a further 600,000 intend to leave the profession by 2027 due to stress, burnout, and retirement. In other words, we are in a vicious circle where understang leads to deterioration, which in turn leads to understang.

In view of these realities and their public health implications, all countries around the world are simultaneously seeking to train, retain and recruit nurses, including, at least for the past twenty years, by engaging in international competition, in particular by recruiting those working in low-income countries (e.g., [START_REF] Brush | Imported Care: Recruiting Foreign Nurses To U[END_REF], [START_REF] Morgan | NHS stang shortages. Why do politicians struggle to give the NHS the sta it needs?[END_REF]), which both the WHO (WHO [2021]) and the International Council of Nurses (ICN [2023]) deplore and try to discourage so as not to further accentuate inequalities between countries and populations.

It is against this backdrop that several countries have begun to consider, without seeing it as a panacea, the creation of ratios of caregivers to patients in hospitals, with the aim of improving both the quality of care and caregivers' working conditions. To our knowledge, state of California (USA) was the rst to enacted a law of this type (Assembly Bill n o 394, Chap. 945, 1999), followed by the state of Victoria (Autralia), Wales (United Kingdom), Scotland (United Kingdom), Ireland, and the state of Queensland (Australia).3 This subject is also currently under discussion in France, thanks to Proposition de loi n o 105 tabled on November 8, 2022 by Senator Bernard Jomier, aimed at modifying the Public Health Code by creating hospital benchmarks in the form of ratios of caregivers per patient.

The aim of our article is to contribute to the analysis of patient/nurse ratios. We propose a way for assessing the number of nurses needed in a hospital department to achieve a given average time for each patient, while complying with legislation on working hours. To this end, we calculate, for a given (potentially random) resource, with random arrivals and exits of patients (which may be of several types), the average time allowed for each patient. To be more precise, let us consider a hospital care unit. During a time step (typically a day), we consider a global care time that medical sta can oer. This global time, which can be random (due to some absenteeism for example) is divided into the patients in the unit. Consider now an additional patient who is coming during the time period under study. The question we raise is to know if the medical sta can deal with this additional patient, given that this patient needs a minimum time of care.

Within the theoretical literature, our paper is in line with the game-theoretic dynamic models with random arrivals that Pierre Bernhard has developed with his co-authors in the context of biology [START_REF] Bernhard | Sharing a resource with randomly arriving foragers[END_REF]) and economics [START_REF] Bernhard | On dynamic games with randomly arriving players[END_REF]][2020][2021]). However, our article diers from these in considering that the resource may be random and that there are dierent types of arrivals (following [START_REF] Biard | Oligopoles avec entrées et types d'entrants aléatoires[END_REF]).

The paper is organized as follows. In Section 2, we present our model in the most simple case where there is only one type of patients (i.e., patients are perfect clones of each other). The arrivals and exits of these patients are random, and we calculate the average time that can be spent on a patient in a hospital service. The results are obtained in closed form for the case of Poisson arrivals.

In Section 3, the previous results are generalized to the case where there are two types of patients, and it helps to understand how to deal with cases where it would be appropriate to consider a greater number of patient types. Finally, in Section 4, we study the inuence of working conditions (summarized as a linear function linking the evolution of the resource from period i to period i + 1) on the average time that can be devoted to a patient.

One prole of patients 2.1 General results

We denote by τ i the available care time during time step i ∈ N in the hospital unit we focus on. The sequence (τ i ) i≥1 can be either deterministic or a sequence of independent random variables with expectation τi := E[τ i ]. In the rst case, we have τi = E[τ i ] = τ i so we can use the common notation τi for both cases. The main question is to compute the care time available for an additional patient who is coming in the unit during the period under study. We assume in this section that all patients are identical, that is to say that the total care time τ i is split equally into the patients. Mathematically, let N i be the number of other patients in this unit at time i ∈ N. So, at time i, there are N i + 1 patients in the unit. Finally, we assume that the (τ i ) i∈N and the (N i ) i∈N are independent, that is to say that the global care time by time step is independent from the number of patients in the unit. The mean care time by patient at time i ∈ N, denoted by ti , is given by

ti = E τ i N i + 1 . ( 1
)
Remark 2.1. A naive approach could state that the mean care time by patient at step i would be equal to the mean global care time divided by the mean number of patients present at time i. However, it is well known by Jensen inequality [START_REF] Jensen | Sur les fonctions convexes et les inégalités entre les valeurs moyennes[END_REF]) that

E[τ i ] E[N i + 1] ≤ E τ i N i + 1 .
As a consequence, the naive approach undervalues the mean care time by patient.

Let us start with N 0 patients in the unit. We assume that N 0 is a random variable over N, independent from all other random variables. At each time i ∈ N * , A i new patients arrive in the unit. We assume that (A i ) i≥1 is a sequence of independent and identically distributed (iid) random variables in N, with common probability-generating function G(x) := E[x A1 ]. At time i ∈ N * , each patient in the unit at time i -1 can leave the unit with probability q. Consequently, given N i-1 , the number of departures at time i ≥ 1, denoted by D i , follows a binomial distribution with parameters N i-1 and q:

D i |N i-1 ∼ Bin(N i-1 , q) .
(2)

Finally, we have:

N i+1 = N i + A i+1 -D i+1 , i ∈ N . ( 3
) Proposition 2.2. Denote by ni := E[N i ] the mean number of patients at time i, for i ∈ N. Let us dene µ := E[A 1 ]
as the mean of the number of arrivals by time step. We have, for all i ∈ N,

ni = p i n0 - µ q + µ q ,
where p = 1q. By consequence,

if n0 = µ q , then ni = µ q for all i ∈ N; if n0 < µ q
, then ni is an increasing sequence;

if n0 > µ q
, then ni is a decreasing sequence.

In the three cases, the sequence is convergent with limit equal to

lim i→∞ ni = µ q .
Proof. From (3), we have

ni+1 = ni + µ -E[D i+1 ] . From (2), E[D i+1 ] = qn i , so ni+1 = (1 -q)n i + µ .
All results come from the fact that (n i ) i∈N is a linear sequence.

This result states that the mean number of patients in the care unit is quite regular and tends to the ratio between the mean number of arrivals at each time over the probability of leaving. The behavior of this mean number only depends on the initial value. If it is lower (resp. higher) than the limit ratio, then the mean number of patients increases (resp. decreases). The particular case where the initial value is equal to the limit gives a constant mean number of patients (see Figure 1 for an illustration).

Figure 1: The mean number of patients is either constant if n0 = µ q , decreasing if n0 > µ q , or increasing if n0 < µ q . See Proposition 2.2.

Remark 2.3. Assume here that N 0 is distributed as A 1 , that is to say that the number of patients at time 0 is distributed as the number of incoming patients during other times. In this case, since n0 = µ < µ/q, the sequence (n i ) i≥0 is increasing.

Lemma 2.4. For i ∈ N, denote by G i (x) := E[x Ni ] the probability-generating function of N i . We have, for all i ∈ N,

G i+1 (x) = G(x)G i (px + q)
Proof. From (2), D i+1 can be rewritten as:

D i+1 = Ni k=1 I k ,
where (I k ) k≥1 is an iid sequence of Bernoulli random variables with parameter q. Thus, from (3),

N i+1 = Ni k=1 (1 -I k ) + A i+1 .
Since A i+1 and N i are independent, we have:

G i+1 (x) = G(x)G i • G 1-I1 (x) ,
where G 1-I1 (x) := E[x 1-I1 ] = px+q since 1-I 1 follows a Bernoulli distribution with parameter 1q = p.

Lemma 2.5. For i ∈ N * , we have

G i (x) = G 0 (p i (x -1) + 1) i-1 k=0 G p k (x -1) + 1 .
Proof. By induction. For i = 0, using Lemma 2.4,

G 1 (x) = G(x)G 0 (px + q) = G 0 (p 1 (x -1) + 1) 0 k=0
G(p 0 (x -1) + 1).

From step i + 1 given step i, we use again Lemma 2.4:

G i+1 (x) = G(x)G i (px + q) , = G(x)G 0 (p i (px + q -1) + 1) i-1 k=0 G p k (px + q -1) + 1 , = G(x)G 0 (p i (px -p) + 1) i-1 k=0 G p k (px -p) + 1 , = G 0 (p i+1 (x -1) + 1)G(x) i-1 k=0 G p k+1 (x -1) + 1 . Since G(x) = G(p 0 (x -1) + 1), and i-1 k=0 G p k+1 (x -1) + 1 = i k=1 G p k (x -1) + 1 ,
we get the result.

Proposition 2.6. We have

t0 = τ0 1 0 G 0 (x)dx,
and

for i ≥ 1, ti = τi 1 0 G 0 (p i (x -1) + 1) i-1 k=0 G p k (x -1) + 1 dx. Proof. From (1), ti = E τ i N i + 1 .
Since τ i and N i are independent,

t i = τi E 1 N i + 1 .
We have

E 1 N i + 1 = E 1 0 x Ni dx = 1 0 E x Ni dx .
Since E[x Ni ] = G i (x) and using Lemma 2.5, we get the result.

As a consequence, to get the mean care time by patient at time i, we only need the probability generating functions of the initial number and the number of arrivals by time step since the integral is easy to get numerically. Note that two dierent distributions for the arrivals give two dierent results even if the mean numbers of arrivals are equal.

The Poisson case

In this subsection, we assume that N 0 is Poisson distributed with parameter λ 0 and the A i 's are Poisson distributed with parameter λ.

Proposition 2.7. We have, for all i ∈ N,

ni = p i λ 0 - λ q + λ q .
By consequence,

if λ 0 = λ q , then ni = λ q for all i ∈ N; if λ 0 < λ q
, then ni is an increasing sequence;

if λ 0 > λ q
, then ni is a decreasing sequence.

In the three cases, the sequence is convergent with limit equal to

lim i→∞ ni = λ q .
Proof. The results are the direct consequence of Proposition 2.2 since the mean number of new arrivals is equal to λ and the mean number of patients present at time 0 is λ 0 .

Proposition 2.8. We have,

t0 = τ0 λ 0 1 -e -λ0 , (4)
and

for i ∈ N * , ti = qτ i (qλ 0 -λ)p i + λ 1 -exp - (qλ 0 -λ) p i + λ q .
(5)

In the case where λ 0 = λ, we have, for all i ∈ N,

ti = qτ i λ (1 -p i+1 ) 1 -exp -λ 1 -p i+1 q . ( 6
)
Proof. Since the A i 's are Poisson distributed, we have -1) .

G(x) = e λ(x
For N 0 ∼ P(λ 0 ), we also have -1) .

G 0 (x) = e λ0(x
So we easily get (4) from Proposition (2.6) and we have

G 0 (p i (x -1) + 1) = e λ0p i (x-1) . Since i-1 k=0 exp λ (p k (x -1) + 1) -1 = exp λ(x -1) 1 -p i q ,
we get

G 0 (p i (x -1) + 1) i-1 k=0 G p k (x -1) + 1 = exp λ 0 p i + λ 1 -p i q (x -1) .
Then, we apply Proposition 2.6 to get (5). Taking λ 0 = λ in ( 4) and ( 5) gives (6).

In the Poisson case, results are totally explicit. In particular, that allows us to study of the behavior of the mean care time by patient over time (see Proposition 2.9). Proposition 2.9. Assume τi is constant over time and denote by τ their common value. We have

if λ 0 = λ q , then ti = qτ λ 1 -e -λ q for all i ∈ N; if λ 0 < λ q
, then ti is a decreasing sequence;

if λ 0 > λ q
, then ti is an increasing sequence.

In the three cases, the sequence is convergent with limit equal to

lim i→∞ ti = qτ λ 1 -e -λ q .
Proof. Letting τ i = τ for all i in Proposition 2.8 gives

ti = qτ (qλ 0 -λ)p i + λ 1 -exp - (qλ 0 -λ) p i + λ q ,
which can been rewritten as ti = τ f (X(i)), with X(i) = (qλ0-λ)p i +λ q > 0 and f (X) = 1-e -X X . Since X ′ (i) = ln(p)p i (qλ0-λ) q and ln(p) < 0, we have that

X(i) is constant if λ 0 = λ q ; X(i) is increasing if λ 0 < λ q ; X(i) is decreasing if λ 0 > λ q . Since f ′ (X) = e -X X-(1-e -X )
X 2 = X+1-e X e X X 2 , the sign of f ′ (X) is the sign of X + 1e X . Since (X + 1e X ) ′ = 1e X which is non positive on (0, +∞), X → X + 1e X is decreasing on (0, +∞), so X + 1e X ≤ 0 for X ∈ (0, +∞) since 0 + 1e 0 = 0. As a consequence, f ′ (X) < 0 for X ∈ (0, +∞) and f is decreasing on (0, +∞), which ends the proof.

In other words, if the mean care time is constant over time, the care time per patient is quite regular and its long-term behavior does not depend on the initial conditions. This initial value determines the monotony of this individual time only through the comparison between the initial mean number of patients and the ratio of the mean new patients per period and the probability of a departure. Figure 2 represents the three behaviors of the sequence ( ti ) i∈N and its convergence.

Figure 2: The mean care time for an additional patient is either constant if λ 0 = λ q , decreasing if λ 0 < λ q , or increasing if λ 0 > λ q . See Proposition 2.9.

Two proles of patients

In this section, only two proles of patients are investigated. Note that the extension to more than two patients types could be obtained using the same methods. Nevertheless, expressions for more than two patients types are quite tedious and that is the reason why it is decided not to present them for the sake of clarity.

General results

We assume throughout this section that two dierent types of patients arrive in the unit. The two proles are called type A and type B. A B-patient needs α more care time than a A-patient. We furthermore assume that α ∈ (1, +∞).

Denote by

N (A) i (resp. N (B) i
) the number of A-patients (resp. B-patients) who are present in the unit at time i ∈ N. We also consider that both

N (A) i i≥1
and

N (B) i i≥1
are iid sequences and independent from one another. As previously, the number of patients are independent from the available care time (τ i ) i∈N . The mean care time for an additional A-patient who arrives at time i ∈ N, denoted by t(A) i , is given by

t(A) i = E τ i (N (A) i + 1) + αN (B) i , (7)
and the time of care for an additional B-patient, denoted by t(B) i , is

t(B) i = E ατ i N (A) i + α(N (B) i + 1
) .

The processes (N

(A) i ) i≥0 and (N (B) i
) i≥0 evolve similarly to that of (N i ) i≥1 in Section 2. Explicitly, let us start with N

G A (A) (x) (resp. G A (A) (x)) the common probability-generation function of the A (A) i 's (resp. A (B)
i 's). Each A-patient (resp. B-patient) present at time i-1 can leave the unit at time i ∈ N * with probability q (A) (resp. q (B) ). Consequently, the number of departures of A-patients (resp. B-patients) at time i, denoted by

D (A) i (resp. D (B) i ), given N (A) i-1 (resp. N (B)
i-1 ), follows a binomial distribution with parameters N (A) i-1 and q (A) (resp. N (B) i-1 and q (B) ):

D (A) i |N i-1 ∼ Bin(N (A) i-1 , q (A) ) , (8)
and

D (B) i |N i-1 ∼ Bin(N (B) i-1 , q (B) ) . (9)
Finally, we have:

N (A) i+1 = N (A) i + A (A) i+1 -D (A) i+1 , i ∈ N , (10)
and

N (B) i+1 = N (B) i + A (B) i+1 -D (B) i+1 , i ∈ N . ( 11
) Proposition 3.1. Denote by n(A) i := E[N (A) i
] the mean number of A-patients at time i and n(B)

i := E[N (B) i
] the mean number of B-patients at time i.

Let ni := n(A) i + n(B)
i the total number of patients at time i. Let us dene µ (A) :=

E[A (A)
1 ] and µ (B) := E[A (B) 1 ] the mean number of arrivals for A-patients and B-patients by time step, respectively. We have

n(A) i = p (A) i n(A) 0 - µ (A) q (A) + µ (A) q (A) ,
and

n(B) i = p (B) i n(B) 0 - µ (B) q (B) + µ (B) q (B) ,
where p (A) = 1q (A) and p (B) = 1q (B) . So

ni = p (A) i n(A) 0 - µ (A) q (A) + p (B) i n(B) 0 - µ (B) q (B) + µ (A) q (A) + µ (B) q (B) .
Letting n goes to innity gives

lim i→∞ n(A) i = µ (A) q (A) , lim i→∞ n(B) i = µ (B)
q (B) , and lim

i→∞ ni = µ (A) q (B) + µ (B) q (B) .
Proof. The proof is similar to the one of Proposition 2.2 since ( 10) and ( 11) are similar to (3).

The behaviors of the mean number of A-patients or B-patients are similar to the mean number of patients in the one-prole case (see Proposition 2.2). In particular, the decrease or the increase of theses numbers only depends on the position of the initial mean numbers with respect to the limit values, which are equal to the ratio between the mean number of arrivals and the probability of leaving for each type of patients. The discussion on the total number of patients is done in Remark 3.2. Remark 3.2. If µ (A) /q (A) < (>)n (A) 0 , then the population of A-patients increases (decreases) with time. The same results hold for the B-patients. Obviously, if the two populations increase (decrease), the total number of patients increases (decreases) too. It is more complicated when one population increases and the other decreases. Indeed, many situations can appear (an increase of the total number of patients all along the time, or a decrease all along the time, or a decrease then an increase or an increase then a decrease). However, we can show (studying i → ni ) that there is one change at most and the monotony changes after time:

T C =       ln -ln(p (B) ) ln(p (A) ) n(B) 0 -µ (B) /q (B) n(A) 0 -µ (A) /q (A) ln p (A) /p (B)       ,
where ⌈x⌉ is the ceiling of x. This conguration can be seen in Figure 3. Proposition 3.3. We have

t(A) 0 = τi 1 0 G (A) 0 (x)G (B) 0 (x α )dx, and t(B) 0 = ατ i 1 0 x α-1 G (A) 0 (x)G (B) 0 (x α )dx. For i ≥ 1, we have t(A) i = τi 1 0 G (A) 0 (p i (x -1) + 1)G (B) 0 (p i (x α -1) + 1) i-1 k=0 G (A) p k (x -1) + 1 G (B) p k (x α -1) + 1 dx, and t(B) i = ατ i 1 0 x α-1 G (A) 0 (p i (x -1) + 1)G (B) 0 (p i (x α -1) + 1) i-1 k=0 G (A) p k (x -1) + 1 G (B) p k (x α -1) + 1 dx.
Proof. From (7),

t(A) i = E τ i (N (A) i + 1) + αN (B) i . Since τ i , N (A) i and N (B) i are independent, t(A) i = τ i E 1 (N (A) i + 1) + αN (B) i .
We have

E 1 (N (A) i + 1) + αN (B) i = E 1 0 x N (A) i +αN (B) i dx = 1 0 E x N (A) i +αN (B) i dx . For t(B) i , from (3.1), t(B) i = E ατ i N (A) i + α(N (B) i + 1) . Since τ i , N (A) i and N (B) i are independent, t(B) i = ατ i E 1 N (A) i + αN (B) i + α .
We have

E 1 N (A) i ) + αN (B) i + α = E 1 0 x N (A) i +αN (B) i +α-1 dx = 1 0 x α-1 E x N (A) i +αN (B) i dx . Since N (A) i and N (B) i are independent, E x N (A) i +αN (B) i = E x N (A) i E x αN (B) i = G (A) (x)G (B) (x α ) .
Since ( 10) and ( 11) are similar to (3), then Lemma 2.5 applies for G (A) (x) and G (B) (x α ), which ends the proof.

Given the probability-generating functions of the number of new incomers of both types and the probability-generating functions of initial numbers of both type, these results give the care time that a new patient of both type could expect at step i, t(A) i and t(B) i . Let us recall here that a B-patient needs twice as much care time as a A-patient. We notice that t(B) i is not equal to α t(A) i but also depends on the dynamics of both numbers of A-patients and B-patients.

The Poisson case

In this subsection, we assume that N Proposition 3.4. We have

n(A) i = p (A) i λ (A) 0 - λ (A) q (A) + λ (A) q (A) ,
and B) .

n(B) i = p (B) i λ (B) 0 - λ (B) q (B) + λ (B) q (
So ni = p (A) i λ (A) 0 - λ (A) q (A) + p (B) i λ (B) 0 - λ (B) q (B) + λ (A) q (A) + λ (B) q (B) .
Letting n goes to innity gives

lim i→∞ n(A) i = λ (A) q (A) , lim i→∞ n(B) i = λ (B)
q (B) , and lim

i→∞ ni = λ (A) q (B) + λ (B) q (B) .
Proof. The results are direct consequences of Proposition 3.1 since the mean numbers of new arrivals are respectively equal to λ (A) and λ (B) and the mean numbers of patients present at time 0 are, respectively, λ for Apatients and B-patients.

Proposition 3.5. We have

t(A) 0 = τi e -λ (A) 0 +λ (B) 0 1 0 exp λ (A) 0 x + λ (B) 0 x α dx, (12)
and

t(B) 0 = ατ i e -λ (A) 0 +λ (B) 0 1 0 x α-1 exp λ (A) 0 x + λ (B) 0 x α dx. (13)
For i ≥ 1, we have

t(A) i = τi e -C (A) i +C (B) i 1 0 exp C (A) i x + C (B) i
x α dx, and

t(B) i = ατ i e -C (A) i +C (B) i 1 0 x α-1 exp C (A) i x + C (B) i x α dx, with C (A) i = λ (A) 0 p (A) i + λ (A) 1 -p (A) i q (A) ,
and Proposition 3.6. Assume τi is constant over time and denote by τ their common value. We have

C (B) i = λ (B) 0 p (B) i + λ (B) 1 -p (B) i q (B) .
lim i→∞ t(A) i = τ e -λ (A) q (A) + λ (B) q (B) 1 0 exp λ (A) q (A) x + λ (B) q (B) x α dx ,
and

lim i→∞ t(B) i = ατ e -λ (A) q (A) + λ (B) q (B) 1 0 x α-1 exp λ (A) q (A) x + λ (B) q (B) x α dx .
In other words, both individual care times are convergent in the case of constant global care time.

Proof. Suppose τi = τ for all i in Proposition 3.5. We have degradation of the working conditions, this time is decreasing afterwards. In Figure 7, the individual mean care time is decreasing and continue to decrease for the three models with a degradation on working conditions but stabilizes when the global mean care time is constant. In both gures, we observe two phases. The rst one ends when the constant case stabilizes and the behavior of the four cases are quite similar. In a second phase, we clearly see dierences.

lim i→∞ C (A) i = λ (A) q (A) , and
lim i→∞ C (B) i = λ (B) q (B) ,
The individual mean care time in the two cases without reaction goes to zero (linearly when the degradation is constant and exponentially when the degradation is brutal). When there exists a reaction, the individual mean care time stabilizes and stay positive. The shape of curves in the second phase is similar to the shape of curves of τi (Figure 5).

Two proles of patients

In this subsection, we consider the setting of Proposition 3.5 where there exist two proles of patients (A and B) and the arrivals are Poisson distributed. Here a B-patient needs twice as much care time as a A-patient. Figure 8 represents the individual mean care time for an additional patient in the case where this patient is of type A or type B. These two individual mean care times are plotted for the four cases of global mean care time. We retrieve the two phases we have mentioned in the previous subsection. We observe the same behaviors as in the one prole case but we can notice that the dierence of individual mean care time between a A-patient and a B-patient reduces with time for the three cases where there exists a degradation of working conditions even in the case where there exists a reaction. Figure 9 gathers the four models for both types.

Conclusion

The main objective of the paper at hand has been to propose a model in order to help assess the number of nurses needed to ensure both healthier caregivers and healthier patients. More exactly, we have proposed a simple model where there is, for everyone, an unknown exact number of patients at each time step. We have assumed that: 1/ all nurses are perfect substitutes (i.e., one hour of a nurse is equal to one hour of another), 2/ patients could be categorized into several types (once and for all), and 3/ entry and exit devices are exogeneous and common knowledge. In this setting we are able to calculate the mean available care time and thereby provide assistance in assessing the potential nursing shortage in a hospital. It should also be noted that we have made the assumption that patient types correspond to pathology types. However, in the event of severe resource scarcity, issues of social choice concerning the distribution of this resource could arise, and patient types could then be reinterpreted as indicating orders of priority. In other words, the model could be reinterpreted in situations where patients need to be sorted.

We left for future research the cases where nurses have dierent skills and must cooperate to care of a given type of patient (i.e., the case where nurses are complement).

To our mind, our analysis has three limits. First, the available care time during time step i (i.e., τ i ) is independent of the number of other patients in the same unit at time i (i.e. N i ), which means that we do not consider situations where the hospital adjusts the number of nurses on the number of patients present at previous periods. Second, our model is only relevant in steady state situations because our random parameters are constant over time. Third, we do not dierentiate between the dierent exit situations (i.e., we do not distinguish the case where the patient leaves because: she is dying, she leaves cured, or she is transported to another service).

Finally, our analysis only focused on the quantity of human resources, but we are aware that there is in fact an interplay between these and, the quality of human ressources (e.g., education level and nursing experience), the material resources, the organization of care services and their cooperation (e.g., [START_REF] Tamata | A systematic review study on the factors aecting shortage of nursing workforce in the hospitals[END_REF]).

  variables in N, independent from all other random variables, and mutually independent. Let A (A) i (resp. A (B) i ) the number of new A-patients (resp. B-patients) at time i. Both sequences (A (A) i ) i≥1 and (A (B) i ) i≥1 are iid sequences of positive random variables in N and mutually independent. Denote by

Figure 3 :

 3 Figure 3: One situation with decreasing mean number of A-patients, increasing mean number of B-patients, and a uctuation in the total mean number (see Remark 3.2).

  distributed with parameter λ(A) , and the A (B) i 's are Poisson distributed with parameter λ(B) .

Figure 4 :

 4 Figure 4: One situation where a B-patient needs twice as much care time as a A-patient (see Proposition 3.5).

Figure 5 :

 5 Figure 5: The four situations for the (τ i ) i∈N during a short period (left-hand side graph) and a long period (right-hand side graph).

Figure 6 :

 6 Figure 6: The four situations for the (τ i ) i∈N during a short period (left-hand side graph) and a long period (right-hand side graph).

Figure 7 :

 7 Figure 7: The four situations for (τ i ) i∈N during a short period (left-hand side graph) and a long period (right-hand side graph).

Figure 8 :

 8 Figure 8: The four situations for (τ i ) i∈N during a short period (left-hand side graph) and a long period (right-hand side graph).

Figure 9 :

 9 Figure 9: The four situations for (τ i ) i∈N during a short period (left-hand side graph) and a long period (right-hand side graph).

There is in particular the reorganization of the hospital with the shift to ambulatory care.

For a positive assessment of this last case, we refer the reader to McHugh et al.[2021].

Proof. The results are the consequences of Propositions 3.3. To get the desired results, we just have to explicit

Since the A (A)

i 's and the A (B) i 's are Poisson distributed, we have

So we easily get ( 12) and ( 13) from Proposition 3.3 and we have

we get -1) . Following the same way gives

which ends the proof.

The Poisson case gives more tractable expressions, which allows us to give the asymptotic behavior (see Proposition 3.6). Figure 4 represents one example of these results where a B-patient needs twice as much care time as a A-patient and the global care time is constant. which ends the proof.

Working conditions impact

In this section, we focus on the quantity τi , that is the available mean care time at step i. We aim to model the impact of a degradation in the working conditions by assuming a decreasing care time along the periods. Mathematically, let τi+1 = aτ i + b , i ∈ N .

Such a sequence is not necessarily decreasing. Let us discuss the cases.

If a = 1, then b has to be negative and greater than -τ 0 /N , where N is the nal time of the study to guarantee the positivity of τi . In this case, we have a linear decrease for the available care time.

If 0 < a < 1, the sequence decreases if τ0 > b/(1a) and converges to b/(1a).

If a > 1, the sequence diverges.

If a < 0, the sequence oscillates.

Only the two rst cases are relevant. The rst case could be seen as a constant degradation without reaction. The second case could model a situation where the degradation is controlled by managers to ensure a minimum service. In this situation, the particular case b = 0 could model an abrupt degradation without reaction (since the limit is zero).

For both cases of one prole and two proles, we are going in this section to compare four dierent situations with respect to the sequence (τ i ) i∈N :

Model 1 A constant situation: τi = τ0 for all i ∈ N.

Model 2 A constant degradation: τi+1 = τi -0.1 for i ∈ N * .

Model 3 A brutal degradation without reaction: τi+1 = 0.99τ i for i ∈ N * .

Model 4 A brutal degradation with reaction: τi+1 = 0.99τ i + 0.2 for i ∈ N * .

The four cases are plotted in Figure 5.

One prole of patients

In this subsection, we consider the setting of Proposition 2.8, that is to say that there exists one prole of patients and the arrivals are Poisson distributed. Let us recall here that in the case where the global care time is constant, Proposition 2.9 states that the mean care time by patient converges and the behavior is regular. In particular, this time is either increasing or decreasing (or constant). For the four models of this section, Figure 6 and Figure 7 depict the two situations for the mean global care time τi . In Figure 6, the individual mean care time is increasing at the beginning but, for the three cases that model a


Marc Deschamps
email: marc.deschamps@univ-fcomte.fr

Pierre Bernhard
email: pierre.bernhard@inria.fr

a discrete time sticky-prices paradox

Keywords: Sticky price, Cournot oligopoly, Dynamic game, Discrete time JEL Classification : C61, C72

come

Introduction

Pricing is clearly at the heart of economic analysis. For a long time analyzed from the viewpoint of price theory, which we now call microeconomics ( [Stigler, 1946]), the desire and need to better microfound macroeconomic models has also led macroeconomists, at least since Keynes, to question this issue ( [Gordon, 1981], [De Vroey, 2016]). As [Mankiw, 1985] and [Romer, 1993] point out, the discussions and oppositions between macroeconomists concerning money, inflation and economic fluctuations are largely centred on the question of price determination by agents.

For some macroeconomists, all nominal prices are perfectly flexible (i.e., all prices correspond to their market equilibrium value), while for others, nominal prices are often sticky (i.e., adjustment is gradual, so there is a difference between the observed price and its theoretical value resulting from market equilibrium). This generally leads the former group to consider that money is neutral in the short term (i.e., there is a dichotomy between the real and monetary spheres), which means that an increase or decrease in the quantity of money has no impact on real economic activity; whereas for the latter group, money is not neutral (i.e., money affects the real economy in the short term) 4 .

Although these debates continue vigorously on the theoretical level, it is noteworthy that today, most macroeconomic models of the dynamic stochastic general equilibrium (DSGE) type are based upon the assumption that firms change their prices only infrequently (see [Woodford, 2003], [Galí, 2015], [Walsh, 2017]). And these models form the academically dominant modeling of the new synthesis ( [START_REF] Goodfriend | The new neoclassical synthesis and the role of monetary policy[END_REF]) and are used daily, as a complement to traditional macroeconomic modeling, by governments, central banks and international agencies ( [START_REF] Christiano | On dsge models[END_REF]).

On the empirical side, beyond older studies (see [START_REF] Klenow | Microeconomie evidence on price-setting[END_REF]] and [Taylor, 1999], for a review of the literature), recent access to numerous and vast microeconomic databases makes it increasingly possible to analyze the issue of price stickiness ( [START_REF] Mackowiak | Implications of microeconomic price data for macroeconomic models[END_REF]). By way of illustration, a recent European Central Bank study of eleven eurozone countries ( [Gauthier et al., 2023]) concluded that: 1/ on average, only 12.3% of prices change each month (8.5% if we exclude sales periods), 2/ differences in terms of price rigidity are limited when comparing countries, and are much greater across sectors, 3/ the median upward price variation over the period 2000-2019 is 9, 6% and 13% downwards (6.7% and 8.7% respectively if sales periods are excluded), and 4/ the distribution of price changes is highly dispersed (14% of price changes are less than 2%, and 10% of price changes are greater than 20%). Generally speaking, it should also be noted that the empirical literature distinguishes between countries with high inflation rates and those with very moderate inflation, the average duration of a price being much shorter in the former than in the latter.

Despite our own limitations in this field, and the fact that it is not our aim here to provide an overview of the issue, we feel that these brief elements are sufficient to pursuade the reader of the theoretical, empirical and political value of studying price stickyness.

The aim of our paper is to contribute to the investigation of the effect of noncontinuous price adjustment in a dynamic Cournot oligopoly, with homogeneous good, n heterogeneous firms and discrete time. We consider an affine inverse demand function, whose parameters (a 0 and b i ) remain constant. There are no demand or supply shocks. We modelize viscosity via the following mechanism 5 : in each period t, a part θ of each firm's output (fixed in time and common to all firms) is sold at the price of the previous period t-1, while the other part of output (1-θ) is sold at the price of period t. This formalization of stickyness seems to us to be one of the simplest imaginable, although exogenous and trivial compared to the literature. In fact, it is like considering that each producer has two warehouses of different sizes (or of the same size), and that at each new period the first warehouse will label the products and sell them at the price of the previous period, while the second warehouse will label the products and sell them at the price of the current period. Alternatively, this assumption could represent the joint time required by each producer to price and notify consumers of the new price, since the period t price applies to both the (1θ)) portion of output in period t and the θ portion of output in period t + 1. Thus, in our model there is a synchronized price adjustment and all firms have the same price duration.

To the best of our knowledge, three articles in the literature are the closest to our own. All three, however, deal with continuous-time dynamic problems where prices evolve continuously. In each case, they analyze the dynamic Cournot-Nash equilibrium, first with open-loop strategies, then with state feedback, the state being the current price. Our paper recovers their common conclusion that the stationary asymptotic price is lower than the Cournot repeated game price but, unlike in ours, in these papers the limit at evanescent stickyness is not the stickinessfree equilibrium. The paper by [START_REF] Fershtman | Dynamic duopolistic competition with sticky prices[END_REF], restricted to the duopoly, offers a study of the asymptotic regime in a framework where dynamics with first-order viscosity makes the current price the state variable, with a quadratic term in the production cost, necessary to avoid a singularity in the continuous time problem. The article by [START_REF] Driskill | Dynamic duopoly with adjustment costs: A differential game approach[END_REF]] is also restricted to the duopoly case and the asymptotic regime, but in their framework what leads to a continuous evolution of prices lies in the fact that the players have for (costly) control the speed of variation of their production rate, and not the rate itself. Thus, their production rate becomes a continuously evolving state variable. The inverse affine law of demand thus produces a continuous variation in price. Finally, the article by [START_REF] Wiszniewska-Matyszkiel | Dynamic oligopoly with sticky prices: off-steady-state analysis[END_REF] considers n players, and is interested in the time trajectory -and not just the asymptotic regime-of prices and production. The dynamics and criteria are the same as in the paper by [START_REF] Fershtman | Dynamic duopolistic competition with sticky prices[END_REF], but they offer a very detailed analysis of trajectories, as well as dependencies in the various parameters.

Our article is organized as follows. In Section 2 we state our dynamic Cournot oligopoly problem with sticky price in discrete time, and explain our assumptions. In section 3 we present the complete solution to our problem, and discuss the dynamics, exhibiting the oscillatory nature of the equilibrium solution and the asymp-totic regime. In section 4 we propose a numerical analysis to compare the sticky price case with the Cournot repeated game. Section 5 analyzes three special cases: the absence of stickiness, monopoly, and Cournot oligopoly when the number n of producers goes to infinity. Section 6 concludes.

The problem 2.1 Cournot dynamic oligopoly with sticky prices

We consider a typical Cournot n firms oligopoly with an affine inverse demand function. Let n be the number of producers, producer i's production be q i , the inverse demand function be characterized by a price a 0 and coefficients b i giving a price P :

P = a 0 - n i=1 b i q i
Each producer i has a linear production cost c i q i . The producers will make an infinite sequence of production decisions q i (t), t ∈ N. But the specificity of this market is that a proportion θ of their production q i (t) will be sold at the previous price P (t -1), while the rest, a proportion (1θ) will be sold at the clearing price P (t) given by the inverse demand function. Therefore, given the appropriate discount factor ρ, player i's profit Π i will be:

Π i = ∞ t=1 ρ t-1 θP (t -1) + (1 -θ) a 0 - n k=1 b k q k (t) -c i q i (t) .
(1)

We seek a Cournot-Nash (dynamic) equilibrium.

The following hypotheses hold on the parameters of the problem :

a 0 > 0 , ∀i , c i < (1 -θ)a 0
so that the a i defined thereafter, are positive. And as in any Cournot model with an affine inverse demand function, we assume that the b i are "sufficiently small" so that realistic productions q i keep P ≥ 0. Furthermore, we will restrict our analysis to the case θ ≤ 1/2. Two reasons lead to this restriction:

1. On the one hand, we have a slightly different interpretation of the same mathematical problem: if the production q i is made at a constant rate q i over the time interval of length 1 between t and t + 1, and the price evolves linearly from x to P during that period, reaching P at a time t + τ < t + 1 and stays there until the end of the period, i.e. time t + 1, then we have the same P (t + 1), and the same profit as expressed by equation ( 1) with θ = τ /2, as in our discrete time problem. Therefore, in this equivalent continuous-time model, θ ≤ 1/2.

2. On the other hand, and more importantly, if θ is too large, the problem may have no solution. To understand this, let us consider the monopoly problem (n = 1) with θ = 1. Then the monopolist may produce a large quantity Q every odd numbered periods (say, years), yielding on even numbered periods a negative6 price P (t -1) which applies for that period when it produces zero, and hence a price a 0 on odd periods. Clearly, its profit will be Π = (a 0c)Q/(1ρ 2 ), hence arbitrarily large. The monopoly problem has therefore no solution in that case. Similar strategies are possible for the n producers model.

Notice also that our formulas will only hold if ρθ 2 /(1θ) 2 < 1, which is ensured (and beyond) by the restriction θ ≤ 1/2.

Notation and preliminary analysis

We will use the following notation:

b i q i = r i , δ := θ 1 -θ , a i := a 0 - c i 1 -θ .
Notice that (1θ)(1 + δ) = 1, so that, e.g., a i = a 0 -(1 + δ)c i . The parameter δ ∈ [0, 1] is an alternative measure of the stickyness, convenient in the calculations, if difficult to interpret in economic terms.

To stress the fact that P (t -1) is the state of the problem at time t, we let P (t -1) = x(t). We may then notice that:

b i 1 -θ Π i = ∞ t=1 ρ t-1 δx(t) + a i - n k=1 r k (t) r i (t) . (2)
and the state dynamics are very simple:

x(t + 1) = a 0 - n k=1 r k (t) . (3)
We have a dynamic game problem with affine dynamics and quadratic payoff. No surprise that we will find a quadratic Isaacs Value function V i (x). We will let

ρW i (x) := ρ b i 1 -θ V i (x) = αx 2 + β i x + γ i . (4)
The fact that α be independent of i will result from the fact that we will succeed in finding such Value functions that satisfy Isaacs' equation.

Further notation used will be:

A = n k=1 a k , D = 1 1 -θ n k=1 c k , and therefore A = na 0 -D , ∆ = n + 1 -2nα , R = n k=1 r k , η i = a i -β i , H = n k=1 η k .

Complete solution

As stated above, we seek an equilibrium with payoffs as in (2) and Value functions as in (4). Isaacs' equation reads as follows:

W i (x) = α ρ x 2 + β i ρ x + γ i ρ = max r i δx + a i - n k=1 r k r i + α a 0 - n k=1 r k 2 + β i a 0 - n k=1 r k + γ i .
(5) This is a concave function of r i . Differentiating and equating to zero, we obtain the equilibrium production r i as:

r i = δx -(1 -2α)R -2αa 0 + η i .
(This expression is still implicit, since R contains r i .) Summing over the i yields, after an elementary calculation:

R = 1 ∆ [n(δx -2αa 0 ) + H] , (6)
and

r i = 1 ∆ [δx -2αa 0 -(1 -2α)H] + η i . (7)
There remains to place this back into equation ( 5) and identify like powers of x.

3.1 Investigation of α and the closed-loop dynamics

Determination of α

Begining with terms in x 2 , we find:

1 ρ α = δ 2 ∆ 2 (n 2 α -2nα + 1) . ( 8
)
We write this equation as:

if δ = 0 , α = 0 , if δ = 0 , 1 ρδ 2 α = f n (α) with f n (α) = n 2 α -2nα + 1 (n + 1 -2nα) 2 .
We observe that 1/ρδ 2 > 1, and furthermore that:

f n (0) = 1 (n + 1) 2 , f n 1 2n = 1 2n , f n 1 2n = 1 , while ∀α ∈ 0, 1 2n + 1 2 , f n (α) > 0 , f n (α) > 0 .
Therefore, in an (y, z) plane, α may be identified as the abscissa of the intersection point of the line z = (1/ρδ 2 )y and the curve z = f n (y). See Figure 1. We know that the slope of the line is larger than one, while the curve z = f n (y) is convex, tangent to the first diagonal at y = 1/2n. There exists therefore one intersection for y < 1/2n, which is the limit of the recursion y(t)/(ρδ 2 ) = f n (y(t + 1)) as t → -∞ starting from y = 0, i.e. the solution that we seek. (Figure 1 easily illustrates the two solutions for y < (n+1)/2n, a third solution is on the decreasing branch of the graph of f n (•) at y > (n + 1)/2n.) A consequence of this graphical representation is that α increases from zero to 1/2n as ρδ 2 increases from zero to one, and that for a given ρδ 2 , α decreases when n increases. Actually, we can even show that 2nα goes to zero as n goes to infinity. (See appendix)

Qualitative behavior of the dynamics

It follows from equations ( 3) and ( 6) that the dynamics under the Cournot-Nash equilibrium strategies are: It follows from the fact that α < 1/2n that ∆ > n. Therefore, nδ/∆ < 1. Hence, for almost all initial conditions, these dynamics oscillate around a long time, asymptotic equilibrium value x:

x(t + 1) = - nδ ∆ x(t) + (n + 1)a 0 -H ∆ . ( 9
) E y 1/2n (n + 1)/2n α T z 1 (n+1) 2 1 2n z = y ρ δ 2 z = f n (y)
x = (n + 1)a 0 -H ∆ + nδ (10)
that we will characterize further later on, when we have calculated H.

Two non-intuitive consequences result from this analysis.

• On the one hand, we insist that this oscillating behaviour is not the result of a trial-and-error process à la Cournot iteration. The actual Cournot-Nash equilibrium srategy yields an oscillation. Our analysis of the extreme case θ = 1 gives an indication of why this may be so.

• On the other hand, although once the prices and productions have reached constant values the stickiness seems to play no role, yet these long term repeated values are not the repetition of the equilibrum values in a game with no stickiness, i.e. x = (a 0 + k c k )/(n + 1).

It may be noticed that, as shown below in subsection 5.2, these somewhat paradoxical facts hold even in the simple case of a one-player game, i.e. a monopoly.

3.2 Investigation of β i and γ i

Coefficients β i and asymptotic price

It is usefull, for more legibility, to introduce yet another short hand notation:

Γ = n 2 -2nα + 1 = ∆ -n(n -1) .
Identifying terms in x in equation ( 4) with the r i as in (7) yields:

β i ρ = δ ∆ 2 (1 + δρ) -2αΓa 0 +2∆(1-nα)a i -2[2nα 2 -(2n+1)α+1]H . ( 11
)
Summing in i, recalling that H = Ai β i , expanding ∆ 2 where it appears without the coefficient δρ, and regrouping terms, we obtain:

∆ 2 δρ + Γ i β i = -2nαa 0 (n 2 -2nα + 1) + 2(n 2 α -2nα + 1)A .
Using again H = Ai β i , and recognizing a term ∆ which appears in the coefficient of A, it follows that

H = [∆ + δρ(n -1)]∆A + 2δρΓnαa 0 ∆ 2 + δρΓ . (12)
Here, the right hand side contains only data of the problem and α that we know how to compute, at least numerically. This is an explicit, although unappealing, formula. It can be placed back into equation ( 11) to get β i and hence η i = a iβ i and place this in equation ( 7) to get the equilibrium strategies. The formulas thus obtained are exceedingly complex and of little interest. If one wants numerical values, the best is to compute H and β i numerically from their respective formulas above.

Asymptotic equilibrium price At this point, we claim that we are able to compute α from formula (8), and we have an explicit formula for H (in terms of α). Therefore, we may compute x with formula (10). We obtain:

x = ∆ 2 (1 + δρ)a 0 + [∆ + (n -1)δρ]∆D ∆ 2 + δρΓ .
See subsection 4 for some numerival values. The dynamics may be written Here a 0 = 10, c = 1, ρ = .95, θ = .5.

x(t + 1) -x = - nδ ∆ (x(t) -x)
5 Particular cases

No stickyness

With no stickyness, we have θ = δ = 0, and consequently, according to our formulas

α = 0 , ∆ = n + 1 , β i = 0 , H = A , η i = a i = a 0 -c i .
We recover the formulas of the classical Cournot-Nash equilibrium with affine inverse demand function. We write them using the shorthand notation n k=1 c k = C and therefore A = na 0 -C, as:

P = a 0 + C n + 1 , r i = a 0 + C n + 1 -c i b i Π i = 1 1 -ρ a 0 + C n + 1 -c i 2 .
In the case where all the production costs coefficients c i are equal, this yields

Π i = 1 b i (1 -ρ) a i n + 1

Monopoly

We now deal with the case n = 1. In that case, α can be calculated exactly. We have ∆ = 2(1α), and equation (8) becomes 4α(1α)ρδ 2 = 0 .

We remember that ρδ 2 is less than one, and let ε = 1ρδ 2 . Then, the smallest root of the above equation yields:

α = 1 2 (1 -ε) , ∆ = 2(1 -α) = 1 + ε , β = ρδ a 1 -(1 -ε)a 0 ρδ + 1 + ε .
The equilibrium dynamics are now

x(t + 1) = - δ 1 + ε x(t) + (1 + ρδ)a 0 + (1 + δ)c 1 ρδ + 1 + ε .
This oscillates around the long term, asymptotic equilibrium

x = (1 + ε) (1 + ε + δ)(1 + ε + ρδ) [(1 + ρδ)a 0 + (1 + δ)c 1 ] .
And finally, the monopoly profit is given by:

Π = 1 1 -ρ 4a 0 [α(1 + δ)c 1 + β] + (a 1 -β) 2 4b(1 -α) .
Numerical values are given in subsection 4 in the column n = 1.

Large number of producers

We may investigate what these formulas say as n goes to infinity. To make things simple, we concentrate on the symmetric case where c i = c for all i, hence A = n(a 0 -(1 + δ)c), and also b i = b. We have already noticed that nα → 0, therefore ∆ ∼ n. It follows that

H ∼ A = n(a 0 -(1 + δ)c) , hence nβ → 0 .
Therefore, R ∼ δx + a and, given that the price x remains bounded, as in the standard Cournot case, all productions go to zero. It is also a simple matter to check that nγ → 0. Hence the producers' profits vanish as well as the cumulative profit of all of them. And finally, the behavior of the asymptotic price is also as in the repeated Cournot case:

x ∼ a 0 + n(1 + δ)c

(1 + δ)n → c .

Conclusion

Our paper presents a simple case of sticky price in a dynamic discrete-time Cournot oligopoly. In this framework, we find the well-known result in the literature that the long-run price is lower than the Cournot repeated game price. On the other hand, in comparison with the continuous case, we establish two new results: 1/ with zero stickiness, the long-run price coincides with the Cournot repeated game price, and 2/ when there is stickiness, the equilibrium solution has an oscillating character. We also show that with a variation in the number of producers (n) or in the share of production sold at the previous price (θ) over time, the optimal solution starts to oscillate again.

It seems clear to us that the discounting of future profits plays an important role in this oscillating character of price trajectories, favoring an increase in profit at time t at the cost of a decrease at time t + 1. Numerical simulations confirm that the discount rate also plays a role in the gap between the asymptotic price and the

Figure 1 :

 1 Figure 1: Determination of α. (Drawing for n = 2, ρδ 2 = 2/3.)

Figure 2 :

 2 Figure 2: Three price trajectories starting from the Cournot price for n = 3, a 0 = 10, c i = 1, and different ρ and θ.

  ρ

Table 1 :

 1 Comparing asymptotic prices and profits for the sticky prices case and repeated Cournot as functions of n.

	.

Table 2 :

 2

	θ	.8	.85	.9	.95	.99
	.48 2.5878 2.5948 2.6002 2.6042 2.6060
	.485 2.5764 2.5825 2.5871 2.5898 2.5903
	.49 2.5644 2.5997 2.5731 2.5742 2.5725
	.495 2.5519 2.5562 2.5581 2.5567 2.5513
	.5	2.5389 2.5417 2.5415 2.5358 2.5196

x for values of ρ and δ close to one with n = 3.

For a broad perspective on theoretical debates on macroeconomics one can refer to[Gaffard, 2018]

We offer in the development another interpretation of the same equations

Or null if we agree that P = max{0, a0i biqi}.

Coefficient γ i and equilibrium profits

To compute equilibrium profits, we still need to evaluate the coefficients γ i . This is obtained by equating terms without x in equation ( 4) with the r i as in ( 7). An explicit expression can be found, but again of little help. As expected, γ i goes to infinity as ρ approaches one.

It finishes to prove that indeed, a Value function of the form (4) can be found that satisfies Isaacs'equation, and therefore that the strategies (7) form a set of a dynamical Cournot-Nash equilibrium strategies.

Starting from a market price P 0 , set x(1) = P 0 , the dynamics (9) gives the sequence of prices under the equilibrium strategies. The total discounted profit of each player is then given by:

The dependence of the final price and profits on θ (or δ) is difficult to assert from these formulas. The following table of numerical values shows that stickyness decreases the price and increases the producers' profits. This is coherent with the rest of the literature on sticky prices.

Some numerical results

We propose here some numerical values aiming to compare the sticky case with the repeated Cournot game. In Table 1, we show on the one hand the asymptotic price P compared to the Cournot price P C , and on the other hand the profits. Since the profit Π i in the sticky case depends on the initial price at time zero, we take it as the Cournot price. The underneath values are for a 0 = 10, c = 1, ρ = .95, and δ = 1, i.e. θ = 1/2:

We give in Figure 2 three price trajectories with the same parameters a 0 and c, n = 3, and different values of ρ and θ.

We also show in Table 2 a nonintuitive phenomenon at very small discount rate. While we expect that the higher the discount rate 1ρ, the more difference we have with the non sticky case, this is not quite so for δ and ρ sufficiently close to one, as the following table shows. We show the asymptotic price P . We have set a 0 = 10 and c = 1 as in Figure 1, and n = 3. We have labeled the columns with ρ and the lines with θ:

Cournot price, a gap that increases with discounting, i.e. as ρ decreases. However, we consider it an open problem that this monotonic growth is no longer true when ρ is very close to 1 and θ close to 1/2.

A Investigation of 2nα

We prove the following theorem:

Theorem 1 For a fixed positive ρδ 2 smaller than one, there exists a unique solution α less than 1/2n of equation ( 8), and 2nα goes to zero as n goes to infinity.

For the sake of clarity, let ω := 2nα, and we call ω the solution sought. Equation ( 8) may be re-written as

Expectedly, we have g n (1) = 1 and g n (1) = 1. We also have g n (0) = 2n/(n+1) 2 . Furthermore, as a simple calculation shows, for all n and all ω ≤ 1, g n (ω) > 0.

As a consequence, g n (•) is a convex function. Hence, ∀ω ∈ [0, 1] , g n (ω) ≤ g n (0) + ω(g n (1)g n (0)) .

Thus, g n (•) being strictly convex,

Let

. so that G n (ω 1 ) = ω 1 /ρδ 2 . Now, observe that ω 1 < 1 so that

The continuous function g n (ω)ω/ρδ 2 is positive for ω = 0 and negative for ω = ω 1 . It follows from the intermediate value theorem that it vanishes at some ω = ω between 0 and ω 1 (and only once because it is convex). We have therefore established that 0 < ω < ω 1 , which goes to zero with g n (0) as n goes to infinity. The intuition for this proof is pictured in the following graphic, which is an enlargment, with a magnification factor 2n, of the lower part of Figure 1.


Simon Gravelle

Roland R Netz
email: rnetz@physik.fu-berlin.de

Lydéric Bocquet
email: lyderic.bocquet@ens.fr

Adsorption kinetics in open nanopores as a source of low frequency noise

Ionic current measurements through solid state nanopores consistently show a power spectral density that scales as 1/f α at low frequency f , with an exponent α ∼ 0.5 -1.5, but strikingly, the physical origin of this behavior remains elusive. Here we perform simulations of particles reversibly adsorbing at the surface of a nanopore, and show that the fluctuations in the number of adsorbed particles exhibit low-frequency pink noise. We furthermore propose theoretical modeling for the time-dependent adsorption of particles on the nanopore surface for various geometries, which predicts a frequency spectrum in very good agreement with the simulation results. Altogether, our results highlight that the low-frequency noise takes its origin in the reversible adsorption of ions at the pore surface combined with the long-lasting excursions of the ions in the reservoirs. The scaling regime of the power spectrum extends down to a cut-off frequency which is far smaller than simple diffusion estimates. Using realistic values for the pore dimensions and the adsorption-desorption kinetics, this predicts the observation of pink-noise for frequencies down to the hertz for a typical solid-state nanopore, in good agreement with experiments.

The fluctuations in the ionic current across a nanopore or nanochannel usually exhibit an anomalous low frequency dependence, in contrast to the expected white noise (flat) spectrum [START_REF] Detcheverry | Thermal fluctuations in nanofluidic transport[END_REF] . Experimentally, the noise spectrum is measured typically to scale as S(f ) ∼ 1/f α at low frequency f , with an exponent α ∼ 0.5 -1.5. This behavior was reported in artificial nanopores and nanotubes with various shapes and materials [START_REF] Siwy | Origin of 1/f α Noise in Membrane Channel Currents[END_REF][START_REF] Chen | Atomic layer deposition to fine-tune the surface properties and diameters of fabricated nanopores[END_REF][START_REF] Tabard-Cossa | Noise analysis and reduction in solid-state nanopores[END_REF][START_REF] Tasserit | Pink noise of ionic conductance through single artificial nanopores revisited[END_REF][START_REF] Smeets | Noise in solid-state nanopores[END_REF][START_REF] Smeets | Low-frequency noise in solid-state nanopores[END_REF][START_REF] Powell | Nonequilibrium 1/f noise in rectifying nanopores[END_REF][START_REF] Powell | Noise properties of rectifying nanopores[END_REF][START_REF] Hoogerheide | Probing surface charge fluctuations with solid-state nanopores[END_REF][START_REF] Wen | Generalized Noise Study of Solid-State Nanopores at Low Frequencies[END_REF] , but also in biological membranes 12,[START_REF] Bezrukov | Examining noise sources at the single-molecule level: 1/f noise of an open maltoporin channel[END_REF] . This low-frequency noise does reduce the signal-to-noise ratio and impedes translocation measurement methods relying on the detection of the modulation of an electrical ionic current accompanying the passage of a macromolecule through a nanopore [START_REF] Howorka | Nanopore analytics: sensing of single molecules[END_REF][START_REF] He | Translocation events in a single-walled carbon nanotube[END_REF][START_REF] Merchant | DNA translocation through graphene nanopores[END_REF] . Such translocation methods are now applied for fast protein detection, genome mapping, or DNA sequencing [START_REF] Shendure | Next-generation DNA sequencing[END_REF][START_REF] Clarke | Continuous base identification for single-molecule nanopore DNA sequencing[END_REF][START_REF] Reisner | Single-molecule denaturation mapping of DNA in nanofluidic channels[END_REF][START_REF] Kowalczyk | Detection of local protein structures along DNA using solid-state nanopores[END_REF][START_REF] Lam | Genome mapping on nanochannel arrays for structural variation analysis and sequence assembly[END_REF][START_REF] Garaj | Molecule-hugging graphene nanopores[END_REF][START_REF] Arjmandi-Tash | Single molecule detection with graphene and other two-dimensional materials: nanopores and beyond[END_REF][24][25] and the low frequency noise reduces the sensitivity of these translocation measurements [START_REF] Chen | Atomic layer deposition to fine-tune the surface properties and diameters of fabricated nanopores[END_REF][START_REF] Smeets | Salt-dependence of ion transport and DNA translocation through solidstate nanopores[END_REF] . Reversly, understanding the origin of this noise is of fundamental importance since it can be a precious source of information [START_REF] Zevenbergen | Stochastic Sensing of Single Molecules in a Nanofluidic Electrochemical Device[END_REF] , for example in electrochemical processes [START_REF] Zevenbergen | Electrochemical correlation spectroscopy in nanofluidic cavities[END_REF][START_REF] Krause | Brownian motion in electrochemical nanodevices[END_REF][START_REF] Singh | Stochastic processes in electrochemistry[END_REF] , on the protonation of ionization sites [START_REF] Bezrukov | Current noise reveals protonation kinetics and number of ionizable sites in an open protein ion channel[END_REF] , the dynamics of adsorption [START_REF] Singh | Stochastic amperometric fluctuations as a probe for dynamic adsorption in nanofluidic electrochemical systems[END_REF] , the transport of sugars [START_REF] Bezrukov | Examining noise sources at the single-molecule level: 1/f noise of an open maltoporin channel[END_REF][START_REF] Nekolla | Noise analysis of ion current through the open and the sugar-induced closed state of the LamB channel of Escherichia coli outer membrane: evaluation of the sugar binding kinetics to the channel interior[END_REF][START_REF] Andersen | Evaluation of the Rate Constants of Sugar Transport Through Maltoporin (LamB) of[END_REF][START_REF] Kullman | Transport of Maltodextrins through Maltoporin: A Single-Channel Study[END_REF] , ATP [START_REF] Rostovtseva | ATP transport through a single mitochondrial channel, VDAC, studied by current fluctuation analysis[END_REF] and antibiotic molecules [START_REF] Nestorovich | Designed to penetrate: Timeresolved interaction of single antibiotic molecules with bacterial pores[END_REF] , the impact of nanobubbles on the ionic current [START_REF] Smeets | Nanobubbles in solidstate nanopores[END_REF] , or as a probe for the measurement of mass flow [START_REF] Mathwig | Electrical cross-correlation spectroscopy: Measuring picoliter-per-minute flows in nanochannels[END_REF] .

The origin of this low-frequency contribution in ionic current has remained elusive. Recent theoretical works have put forward the role of ion correlations [START_REF] Zorkot | The Power Spectrum of Ionic Nanopore Currents: The Role of Ion Correlations[END_REF] as well as hydrodynamic interactions [START_REF] Zorkot | Current fluctuations in nanopores: The effects of electrostatic and hydrodynamic interactions[END_REF] , as a source of low frequency noise. On the experimental side, low frequency noise measured in carbon and boron-nitride nanotubes 25 showed dependence on the surface charge. This suggests that surface effects and charge adsorption on the channel surface play an key role for the low-frequency spectrum. A usual result of the experiments is that the level of current noise obeys the so-called Hooge's low S(f ) = A × I 2 /f α , with I the mean current, and the coefficient A being inversely proportional to the number of charge carriers in the sample [START_REF] Smeets | Noise in solid-state nanopores[END_REF]25,[START_REF] Hooge | 1/F Noise Is No Surface Effect[END_REF] . This points to the fluctuations of the channel conductance, G, as the source of the noise spectrum. Indeed the current fluctuation δI is related to the conductance fluctuations δG according to δI 2 (f ) = δG 2 (f )

G 2
× I 2 , in line with Hooge's law [START_REF] Hoogerheide | Probing surface charge fluctuations with solid-state nanopores[END_REF]25 . Now, the ionic conductance is proportional to the number of charge carriers in the pore, so that its fluctuations originate in the fluctuations of the number of 'free' ions in the channel.

In this work we accordingly explore the role of adsorption-desorption kinetics of particles (charged and neutral) on the surface of nanopores: individual particles diffuse in the bulk of a nanopore, and possibly in the reservoirs, and experience adsorption-desorption processes on the nanopore surface. Using molecular simulations and theoretical modeling, we calculate the frequency spectrum of free and bounded particles in the nanopore. We explore how the geometry affects the low frequency spectrum, as well as the impact of the presence of reservoirs connected to the nanopore.

Simulations.

We start with the description of the process using Molecular Dynamics (MD) simulations. MD simulations are performed using LAMMPS [START_REF] Plimpton | Fast Parallel Algorithms For Shortrange Molecular-dynamics[END_REF] . A typical system is shown in Fig. 1 A and consists of a cylindrical nanopore (here with radius R n = 2 nm and length L n = R n ) connected to large half-spherical reservoirs (here with radius H n = 20 nm). The index 'n' stands for 'nanopore'. Charges are homogeneously distributed among the wall atoms at the surface with a surface charge Σ. A number N = 2πR n L n ×Σ/e of monovalent counterions is initially placed in the nanopore in order to assure electroneutrality (here N = 5 corresponding to Σ 32 mC/m 2 ); see Fig. 1 B. Monovalent ions of concentration c s = 15 mM are added to the reservoir. The temperature of the ions is fixed to 300 K thanks to a Langevin thermostat with a damping time τ = 23 fs 45 , leading to a diffusion coefficient D ≈ 1.5 × 10 -5 cm 2 /s. The equations of motion are integrated using the velocity Verlet algorithm [START_REF] Martys | Velocity Verlet algorithm for dissipative-particle-dynamics-based models of suspensions[END_REF] , with a time-step equal to 2 fs. Long-range Coulombic in-teractions are calculated by means of a particle-particle particle-mesh (PPPM) solver [START_REF] Darden | Particle mesh Ewald: An N log(N) method for Ewald sums in large systems[END_REF] . In addition, particles interact through a Weeks-Chandler-Andersen potential (σ = 3 Å), and the reversible adsorption of particles at the nanopore surface is introduced by adding a supplementary harmonic potential between fluid and wall atoms. A counter-ion forms a bond with a wall atom when the distance of their respective centers is lower than σ. The particle is then trapped to the surface site by a harmonic potential; the bond however breaks when the distance between the particle and the wall atom's center is higher than σ. The particle may thus escape thanks to thermal fluctuations. The depth of the harmonic potential, similar to the adsorption energy, is equal to = 0.5 Kcal/mol. This leads to an adsorption duration T distributed as P (T ) = λ exp(-λT ), with a desorption frequency λ ≈ 0.8 GHz. More details are given in the Appendix section. We record the number of bounded counter-ions N B at the cylinder surface as a function of time, as well as the total number of free (nonbounded) ions within the nanopore, N F . Let us mention that the value for λ in the simulations was chosen to be substantially larger than its experimental counterpart, typically in the range λ exp ∼ 10 4 Hz 31 . This allows us to observe a sufficient number of desorption events within the range of the simulation time for proper statistics and quantify the influence of the adsorption-desorption kinetics in comparison to the theoretical predictions.

We calculate the Power Spectral Density (PSD) of both signals N B (t) and N F (t) [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF] . As shown in Fig. 1 C, the PSD of the number of bounded particles N B scales as 1/f α with α ≈ 1.5, for frequency f in between 5 • 10 -3 and 5 GHz (Fig. 1 C). Note that the fluctuations of δN B are proportional to the fluctuations of the surface charge density δΣ. Fig. 1 D then shows the PSD of the number of free particles N F in the pore. The latter exhibits a frequency dependence behaving as 1/f α with α ∼ 0.7 in between 5 • 10 -3 and 0.5 GHz, and α ∼ 1.2 for larger frequencies.

We compare these results with those obtained by simulations of three simpler model systems: (i) We first considered a neutral system, for which the ions are replaced by a number N = 100 of neutral particles and a charge-free surface. As shown in Fig. 1 C, (open symbols), this leads to almost identical PSD for the number of bounded particles. However, we find a substantial difference on the PSD of the number of free particles as compared to the charged case, since the scaling of the PSD for the neutral system exhibits a slightly lower exponent (α ∼ 0.5) as compared to the charged system (Fig. 1 D). (ii) We also explored simulations of non-interacting particles, since correlations between ions were shown to be a potential source of pink noise [START_REF] Zorkot | The Power Spectrum of Ionic Nanopore Currents: The Role of Ion Correlations[END_REF] . We used a conventional Brownian dynamics algorithm with overdamped dynamics, implemented with a homemade C++ code (Method section). Strikingly, Brownian simulations lead to results very similar to neutral MD simulations (Fig. 1 C D), indicating that the observed pink noise is not a consequence of the correlations between the particles induced by the WCA potential. In addition, because Brownian simulations allow for the exploration of long duration at a low computational cost, the results reveal a plateau of the PSDs at frequencies below f 0 = 4 • 10 -3 GHz, for both PSDs of the bounded N B (t) and free N F (t) particles. We will refer to f 0 as the low frequency cut-off. (iii) Finally, we explored the PSD in an infinite cylinder, with a radius equal to that of the nanopore (R c = 2 nm, where 'c' stands for 'cylinder'), periodic boundary conditions, and an identical adsorption energy (Fig. 2 A). Note that in the infinite system, δN F = -δN B and their respective PSDs are identical. For high frequencies the PSD of the bounded particle number show the same scaling as the one obtained with the nanopore connected to reservoirs, i.e. 1/f α with α ≈ 1.5. However it exhibits a plateau for frequencies below 0.1 GHz (Fig. 1 C, pale blue). The low frequency noise is therefore strongly affected by the presence of the reservoirs, a behavior which we will analyze more in details below.

Modeling of the adsorption kinetics.

In order to isolate the origin of the pink noise detected in both MD and Brownian simulations, let us focus first on neutral particles reversibly adsorbing at the surface of an effectively infinite cylinder (Fig. 2 A). Our approach extends on the theoretical developments on number fluctuations in a channel by Bezrukov et al. [START_REF] Bezrukov | Particle number fluctuations in a membrane channel[END_REF] , as well as on the analysis by Levitz et al. on random flights on planar surface [START_REF] Levitz | Random flights in confining interfacial systems[END_REF][START_REF] Levitz | Molecular intermittent dynamics of interfacial water: probing adsorption and bulk confinement[END_REF] . We define a state variable n i B (t) which is equal to 1 when a particle is bounded on the surface, and 0 otherwise. The number of bounded particles at any time t follows as N B (t) = N i=1 n i B (t). We assume that particles are uncorrelated among each other. The auto-correlation function of the total number of adsorbed particles can be written as a sum

N B (0)N B (t) = N 2 p 2 B + N p B C BB (t) -N p 2 B , (1)
where p B is the probability for a particle to be adsorbed and C BB (t) the single-ion correlation function, i.e. the probability for an ion to be adsorbed at t given that it is adsorbed at t = 0. The correlation function C BB can be calculated exactly in terms of the successive adsorption-desorption-diffusion processes of the particles. Explicitly this takes the form in the present case:

C BB (t) = ∞ n=0 ∞ 0 dt e Q(t e ) n m=1 ∞ 0 dt m P (t m ) ∞ 0 dt m J c (t m ) × δ t e + n k=1 (t k + t k ) -t , (2)
where the outer sum counts the number of desorption and re-adsorption events, P (t) is the distribution of the desorption duration, P (t) = λ exp(-λt), and Q(t) is the survival distribution corresponding to the probability for a particle adsorbed in t = 0 to remain adsorbed at a time t, given by

Q(t) = ∞ t P (t )dt . (3)
Finally, J c (t) is the first-return distribution, defined as the probability for a particle desorbing at time t = 0 to return to the cylinder surface at t. The Laplace transform of C BB (t) allows to perform the sum as a geometric series, leading to the simple result

CBB (ω) = Q(ω) 1 -P (ω) Jc (ω) , (4)
where P (ω) = λ (λ + ω) -1 and Q(ω) = (λ + ω) -1 are respectively the Laplace transforms of P (t) and Q(t), and ω = 2πf .

For the cylinder geometry, the first-return distribution Jc (ω) is calculated by solving the diffusion equation in cylindrical coordinates

(∂ t -D∆ r )G c (r, t|r 0 ) = 0, (5)
where G c is the radially symmetric Green function and

∆ r = r -1 ∂ r r∂ r ; Jc (ω) is related to G c via the relation Jc (ω) = -D∂ r Gc (R n , ω|r 0 = R n ).
The diffusion equation is complemented by the initial condition G c (r, t = 0|r 0 ) = δ(r -r 0 )/2πr, and the boundary condition on the pore surface which writes

-D∂ r G c (r, t|r 0 )| r=Rc = k G c (R c , t|r 0 )
, where we introduced the effective adsorption rate k (m/s). This leads to an expression for the first-return distribution for an infinite tube of radius R c Jc (ω

) -1 = 1 + √ ωD k I 1 (R c ω/D) I 0 (R c ω/D) , (6)
where I 0 (x) is the modified Bessel functions and 52 . Further details of the derivation are given in the Supplementary Information. The PSD follows as

∂ x I 0 (x) = I 1 (x)
S(ω) = CBB (iω) + CBB (-iω), (7)
with i the imaginary unit, and CBB from Eq. ( 4). We compare this prediction to the results of Brownian simulation in Fig. 2 D, using three values of the effective desorption frequency, λ = 0.01, 10 and 100 GHz and a radius R c = 20 nm. Note that the value for the effective adsorption rate k = 22 m/s is directly measured from the simulations: k is obtained in terms of the average surface concentration, c s (mol m -2 ), and bulk concentration in particles, c b (mol m -3 ), according to k = λc s /c b . Overall, the agreement is shown to be excellent. For a large desorption frequency, λ = 100 GHz, corresponding to particles quickly returning to diffusing after being adsorbed at the surface, the PSD shows three distinct regimes: a plateau at low frequency (f < 4 • 10 -3 GHz), a decrease in 1/f α with α ≈ 1.5 at high frequency (f > 10 GHz), and in between a decrease in 1/f α with α ≈ 0.5. As the desorption frequency λ decreases, the intermediate regime eventually disappears and only two distinct scaling regimes remain (Fig. 2 D). Note that a similar impact of the desorption frequency on the scaling of the PSD of confined molecules alternating between diffusion and adsorption has already been reported by Levitz et al. [START_REF] Levitz | Molecular intermittent dynamics of interfacial water: probing adsorption and bulk confinement[END_REF] . Regardless of the desorption frequency λ, the low frequency cut-off f 0 is given by the diffusive characteristic frequency f 0 = f D ≈ D/R 2 c , as indicated by both Brownian and neutral MD simulations (Fig. 2 E).

The results obtained for cylinders can be generalized to alternative geometries (Fig. 2 B,C), with planar and spherical confinement and adsorbing surfaces. On the theoretical side, for the sphere, the previous framework still holds, with the first-return distribution Jc (ω) replaced by those calculated for spherical geometry Js (ω) (Supplementary Information, Eq. (S45)). For the planar geometry, the equations are slightly more complicated because one must account for particles desorbing from one surface and re-adsorbing to the other (Supplementary Information, Eq. (S59) to Eq. (S86)). As before, the agreement with the simulations is found to be excellent for both the planar geometry of Fig. 2 B, with the particles diffusing in between two adsorbing plates, see Fig. 2 F, and for the spherical geometry of Fig. 2 C, see Fig. 2 G.

The study of diffusion of particles within a confined geometry with adsorbing boundaries, either in planar, cylindrical, and spherical, all lead to the same results: the fluctuations of the number of bounded particles exhibits pink noise, with a PSD scaling as 1/f α with α ∈ [0.5, 1.5], for frequencies above a cut-off fixed by the inverse duration of diffusion trajectories

f 0 = f D ≈ D/R 2
x , where R x is the characteristic dimension of the considered system, i.e. cylinder or sphere radius, or plate separation. The origin of the pink noise is the diffusive return trajectories to the surfaces, combined with the reversible adsorptions. However, the value of the cut-off frequency remains extremely large: indeed even for a large nanopore with radius R n = 100 nm, f 0 = f D ∼ 10 4 Hz for typical ions. This is in contradiction with experimental results, which exhibit pink noise down to extremely small frequencies in the range of Hz. Furthermore, our previous numerical results for the complete system including reservoirs exhibited pink noise down to frequencies which are far smaller than f D , see Fig. 1 C D. We thus extend the analysis to explore the effect of reservoirs connecting the pore.

Reservoir effects.

In order to obtain insight on the cut-off frequency in the presence of reservoirs, we study a simplified geometry where an adsorbing system is connected to a single large reservoir. Accordingly we focus on a system made of two concentric spheres of respective radius R cs and H cs , with H cs > R cs , as sketched in Fig. 3 A. The surface of the inner sphere is adsorbing, and the surface of the outer sphere is reflecting. Particles are located in between R cs and H cs . The limit R cs H cs is expected to mimic, at least roughly, a nanopore connected to large reservoirs. This geometry captures the elementary process at play when ions escape the nanopores, while keeping analytical calculations tractable.

In order to calculate the analytical expression for the PSD, we solve the diffusion equation in spherical coordinates, (∂ t -D∆ r )G s (r, t|r 0 ) = 0, where G s is the radially symmetric Green function and ∆ r = r -2 ∂ r r 2 ∂ r . The initial condition G s (r, t = 0|r 0 ) = δ(r -r 0 )/4πr 2 and the boundary condition on the surface at r = R cs now writes (k + D∂ r )G s (r, t|r 0 )| r=Rcs = 0. In the limit of large reservoir, H cs R cs , we obtain the first return distribution of particles onto the inner sphere

Jcs (ω) -1 = 1 + D kR cs + √ ωD k . ( 8
)
Further details are given in the Supplementary Information. The calculation of the PSD of the number of bounded particles at the inner sphere surface then follows along the same lines as above. The Laplace transform of the correlation function writes

CBB (ω) = Q(ω) 1 -P (ω) Jcs (ω) , (9)
with Jcs from Eq. ( 8), and the PSD S(ω) is deduced as S(ω) = CBB (iω) + CBB (-iω). To assess the validity of the calculations, we compare to numerical simulations in Fig. 3. The inner sphere radius is R cs = 8 nm, and the outer sphere radius is H cs = 10 × R cs , which is large enough to ensure that the outer sphere has no impact on the measurement, as confirmed by convergence tests. We performed simulations of uncorrelated Brownian particles. For a large desorption frequency, λ = 100 GHz, the PSD of the number of bounded particles at the surface of the inner sphere, N B (t), is qualitatively similar to the PSD obtained with confined geometry. The PSD possesses three distinct regimes: a flat PSD at low frequency, a PSD ∼ f -3/2 at high frequency, and an intermediate regime with PSD ∼ f -1/2 (Fig. 3 B). As the desorption frequency decreases, the intermediate regime disappears. The agreement of the theoretical predictions with the simulation results is excellent, see Fig. 3 B.

But these results are mostly interesting to obtain further insights on the cut-off frequency, which is the question at stake here. As highlighted by the simulations, see Fig. 3 B, a crucial difference between the concentric spheres and confined system is that the low frequency cut-off is pushed down to much lower frequency as compared to the diffusive one when the desorption frequency decreases. The diffusive cut-off f D = D/(2πR 2 cs ) takes the value f D ≈ 4 × 10 -3 GHz, but for the lowest desorption frequency λ, the scaling behavior PSD ∝ f α is preserved down to ∼ 10 -5 GHz (Fig. 3 B). This echoes the results obtained with a nanopore connected to large reservoirs in Fig. 1 C D. Physically, our results indicate that the typical duration of the correlations occurring in the system is much longer in a system connected to reservoirs, as compared with a confined geometry.

The dependency of the cut-off frequency ω 0 as a function of the desorption frequency λ, 'pore size' R, and effective adsorption coefficient k is shown in Fig. 3 C. This is compared to the analytical prediction obtained from the calculation of the PSD above. The predicted low frequency cut-off f 0 = ω 0 /2π is deduced from

f -1 0 = f -1 D + f -1 λ , (10)
where f D = D/(2πR 2 cs ) and

f λ = λ/2π (1 + kR cs /D) , (11)
where the index 'λ' refers to the dependency in the desorption frequency λ. These predictions do match the simulation results with varying desorption frequency λ and varying inner sphere radius R cs (Fig. 3 Ca, Cb). In order to test further the validity of Eq. ( 10), we also explored the effect of an artificial increase of the effective desorption rate k = λc s /c b . To do so, we have added a radial force to the particles forcing them to return to the surface of the inner sphere, hence increasing the average fraction of particles adsorbed at the inner sphere surface (Appendix section). For a variation of k from 1.5 m/s to 980 m/s, the agreement between simulation and Eq. ( 10) is also excellent.

Back to the nanopore. These results allow us to discuss more in depth the results for the realistic nanopore in Fig. 1. Indeed, while the above geometry is simplistic, it does capture the elementary process of particle return which is at play in the reservoirs of the realistic nanopore in Fig. 1. The above Eq. ( 10) predicts a cut-off of f 0 = 4 × 10 -3 GHz for the system in Fig. 1, which agrees well with the measured cut-off. Note however that in Fig. 1, the cut-off is only clearly visible for the Brownian simulations and not from MD due to numerical limitations at low frequencies.

To validate further Eq. ( 10) when applied to the case of nanopore, we performed MD simulations using a charged nanopore of varying radius R n = L n , with H n = 10R n . The surface charge is Σ = 60 mC/m 2 and salt concentration is c s = 20 mM. In order to make the cut-off f 0 detectable given the limitation of the simulation, we use a high desorption frequency λ ≈ 100 Ghz. The cut-off predicted by Eq. 10 gives an excellent agreement with simulations, using the value of k = 22 m/s (Fig. 4 A). MD simulations of neutral systems and Brownian simulations of uncorrelated particles also fall on the same curve. To further confirm the validity of Eq. 10, we also vary the desorption frequency λ for a nanopore of radius R n = 1 nm and find again an excellent agreement between simulations and Eq. ( 10) (Fig. 4 B).

Finally we report the exponent of the PSD for both the number of bounded N B and free N F particles, where the exponent α is defined from PSD ∼ f -α . MD results show that α varies strongly with the bulk salt concentration c s , in particular the exponent of the PSD of N F ranges from α ∼ 0.15 to ∼ 1.3 (Fig. 4 C). The exponent also depends on the desorption frequency λ (Fig. 4 D), which echoes the results obtained with the two concentric spheres (Fig. 3 B). Finally, we did vary the surface charge from 60 to 360 mC/m 2 , but did not observe significant variation of the exponent α, nor significant variation of the cut-off f 0 (Supplementary Information, Fig. S1). However, the measured increase of k was only by a factor 2 in between Σ = 60 and 360 mC/m 2 , which is expected to impact f 0 by about 10 % according to Eq. ( 10), hence smaller than the error bars. Note that we did however observe an increase of the level of noise with increasing surface charge Σ (Supplementary Information, Fig. S1), in qualitative agreement with the experimental observations in Ref. 25.

Discussion.

In this work we have shown that the fluctuations associated with particles reversibly adsorbing at a solid surface lead to a signal with a PSD scaling as 1/f α , with an exponent α typically in the range 0.5 to 1.5. We recall that the conductance fluctuations are expected to be directly proportional to the PSD of the free ions the fluctuations in the number of free ions, so that the results for the free particles PSD apply to the conductance as well. A crucial remark though is that this low-frequency scaling behavior is observed only above a minimum frequency cut-off f 0 . Our predictions show that f 0 depends not only on the system size and particles diffusion coefficient, but also on the desorption energy and adsorption rate. When the adsorption-desorption process is limiting (k D/R c ), the cut-off frequency predicted from Eq. ( 10) is typically

f 0 = ω 0 2π ≈ λD 2πkR c . (12)
Based on these results, we are now in position to discuss the experimental observations. While the prediction of 1/f α noise is indeed in agreement with the experiments, it is essential to assess that the range of observation for also matches the experimental ones whether the predicted range for the pink noise also matches experimental ones. In terms of orders of magnitude, simulations showed that a cut-off frequency of about 0.004 GHz is obtained for λ ≈ 0.8 GHz. While this is still very high compared to the experimental range, one should now estimate the cutoff frequency f 0 for more realistic values of the adsorption-desorption kinetics. First, the reported values for the experimental desorption frequencies λ are rather in the range λ = 10 4 -10 5 Hz, as measured for the protonization of a transmembrane ionic channels [START_REF] Bezrukov | Current noise reveals protonation kinetics and number of ionizable sites in an open protein ion channel[END_REF] , as well as on the surface of solid-state SiN nanopores [START_REF] Hoogerheide | Probing surface charge fluctuations with solid-state nanopores[END_REF] . Also based on the association rate constant k R ∼ 10 11 -10 12 M -1 s -1 infered for these experiments, one may obtain an estimate for the parameter k (in m.s -1 ) as k = k R × Γ, with Γ the number of reactive sites per unit surface (in m -2 ). Using the values for SiN 10 (Γ = 0.4 nm -2 ), one gets k ∼ 100 m/s. Altogether, for a typical pore of size R n ∼ 10 nm, Eq.( 12) predicts a cut-off frequency as low as 3 Hz for a solid state nanopore. Similar estimates for nanotubes such as those studied in Ref. 25 lead to an even smaller cutoff in the range of a fraction of Hertz. These estimates show that for the experimentally relevant values of the surface adsorption-desorption coefficients, the cut-off frequency may reach sub-Hertz values and thus lower than the range of frequencies investigated experimentally. In other words the low frequency scaling for the ionic noise S(f ) ∼ f -α is expected to be observed down to very low frequencies in solid-state nanopore and no white noise regime is observed. However, it is interesting to note that the white-noise regime was actually observed in some measurements performed on protein ion channel [START_REF] Bezrukov | Current noise reveals protonation kinetics and number of ionizable sites in an open protein ion channel[END_REF] , with a pink noise behavior which is observed only above a cut-frequency of about 5 kHz. In this system, pink noise is also expected to originate from the reversible adsorption of protons at ionizable sites. Using the reported values for λ and k R for this system, one may deduce that k ∼ 10 m/s. Using a typical radius of R n ∼ 0.5 nm for these pores [START_REF] Menestrina | Ionic channels formed by Staphylococcus aureus alpha-toxin: Voltage-dependent inhibition by divalent and trivalent cations[END_REF] , one gets f 0 ∼ 1.5 kHz from Eq. ( 12): this is indeed in the same order of magnitude as the experimentally observed cut-off. These experimental results hence support the importance of surface adsorptiondesorption processes as the origin of low-frequency pink noise. As a further confirmation of surface effects, the amplitude of the noise in carbon nanotubes was measured in Ref. 25 to be strongly dependent on the surface charge. This result is also in agreement with results from the simulation, see supplementary Fig. S1 Altogether our theoretical description demonstrates that the counter-intuitive noise spectrum observed in nanopores takes its origin in the adsorption-desorption process on the pore surface coupled in a crucial way to the long-lasting excursions of the ions into the reservoirs.

(Fig 2 A), or two concentric spheres (Fig 3 A). The particles do not interact with each other (particle interaction potentials set equal to 0). The simulation time is iterated with constant time step ∆t = 1 ps, and particles positions are iterated following Brownian dynamics (also called overdamped Langevin dynamics) procedure: a Gaussian random number r is generated, and the particle jumps are calculated as ∆r = √ 6D∆t r, with D = 1.5 • 10 -5 cm 2 /s the particles diffusion coefficient. When a particle reaches an adsorbing wall, it adsorbs with a probability p abs = 1 for a duration T distributed as P (T ) = λ exp(-λT ) with λ the desorption frequency. An adsorbed particle is maintained fix at the surface during T . Data acquisition starts after an equilibrium phase. In the case of concentric spheres, if specified, an additional radial force is added to the particles. The additional force is included by adding a deterministic radial velocity to the particles.

Figure 1 :

 1 Figure 1: (A) Sketch of a cylindrical nanopore of length Ln, radius Rn, and connected to two half spherical reservoirs of radius Hn. The disks represent co-ions (green) and counterions (blue). (B) Snapshots of the molecular dynamics system 40 . The blue and green disks are the ions and the red disks are the traps. Solid flat walls are highlighted in gray. (C) PSDs of the number of bounded particles at the adsorbing surface as function of the frequency f for a charged system (full disks), as well as for a neutral system (open squares), uncorrelated Brownian particles (black dots), and effectively infinite cylinder (full squares). The dashed line highlights the exponent as f -α . PSDs have been vertically shifted for clarity. (D) Corresponding PSDs of the number of free particles within the nanopore. The PSD for the infinite cylinder is not represented.

Figure 2 : 1 0

 21 Figure 2: (A) Particles (blue disks) within an effectively infinite cylinder of radius Rc and length Lc with adsorbing boundary (red surface). (B) Particles in between two effectively infinite flat planes of lateral dimension Lp separated by a distance 2Rp. (C) Particles within a sphere of radius Rs. (D) PSD of the number of bounded particles at the adsorbing surface of an infinite cylinder as function of the frequency f for a tube radius Rc = 20 nm and a desorption frequency, respectively from bottom to top, λ = 100 (red disks), λ = 10 (gray triangles) and λ = 0.01 GHz (blue squares). Dashed lines are guides to the eyes ∝ f -α , full lines are the analytical calculations (Eq. (7)), and dotted line is f0, see text for details. PSDs have been vertically shifted for clarity. (E) Measured inverse cut-off frequency f -1 0 as function of the cylinder radius Rc for three different values of λ, the black line is R 2 c /D with D the diffusion coefficient of particles. Green stars are neutral MD simulations. (F) PSD for a planar system with separation 2 × Rp = 40 nm and a desorption frequency λ = 100 GHz. (G) PSD for a sphere radius Rs = 20 nm and a desorption frequency λ = 100 GHz.

Figure 3 : 1 0

 31 Figure 3: (A) Concentric spheres of respective radius Rcs and Hcs, with adsorbing boundary in red and reflecting boundary in gray. (B) PSD of the number of bounded particles at the adsorbing surface of concentric spheres as function of the frequency f for spheres of respective radius Rcs = 8 nm and Hcs = 80 nm, and a desorption frequency, respectively from bottom to top, λ = 100 (red disks), λ = 1 (gray triangles) and λ = 0.01 GHz (blue squares). Dashed lines are guides to the eyes PSD ∝ f -α , full lines are the analytical calculations (Eq. (9) with the result (8)). PSDs have been vertically shifted for clarity. (Ca) Inverse of the low frequency cut-off f -1 0 as measured from simulations for Rcs = 8 nm (symbols), as a function of the inverse desorption frequency λ -1 . The full line is Eq. (10), with k = 22 m/s. (Cb) Inverse of the low frequency cut-off f0 as function of the inner sphere radius Rcs for λ = 100 (red disks), λ = 1 (gray triangles) and λ = 0.01 GHz (blue squares), and for k = 22 m/s. (Cc) Inverse of the low frequency cut-off f0 as function of the effective adsorption boundary k for Rcs = 8 nm and λ = 1 GHz.

Figure 4 :

 4 Figure 4: (A) Inverse of the low frequency cut-off f0 as function of the nanopore radius Rn, for a charged nanopore with surface charge Σ = 60 mC/m 2 , salt concentration cs = 20 mM, and desorption frequency λ = 100 GHz (full symbols). We also report results obtained with a neutral nanopore (open symbols), and using Brownian simulations (stars). The solid line is the prediction in Eq. (10) with k = 22 m/s. (B) Inverse of the low frequency cut-off f0 as a function of the inverse desorption frequency λ -1 for a nanopore with radius Rn = 1 nm. (C) Slope α of the PSD of the number of free particle (blue disks) and number of bounded particle (red squares), for a charged nanopore with surface charge Σ = 60 mC/m 2 , desorption frequency λ = 0.8 GHz, and nanopore radius Rn = 2 nm. Lines are guides to the eye. (D) Slope α as a function of the desorption frequency for a charged nanopore with surface charge Σ = 60 mC/m 2 , nanopore radius Rn = 1 nm, and salt concentration cs = 20 mM. Lines are guides to the eye.

Acknowledgments

LB and RRN acknowledge support from ANR grant Neptune and DFG grant NE 810/11. LB acknowledges support from the European Union's Horizon 2020 Framework Program/European Research Council Advanced Grant 785911 Shadoks. RRN acknowledges support from DFG grant SFB 1078. The authors / SG thank(s) Laurent Joly and Francois Detchevery for valuable discussions.

Appendix

Molecular dynamics simulations

We performed Molecular Dynamics (MD) simulations using LAMMPS [START_REF] Plimpton | Fast Parallel Algorithms For Shortrange Molecular-dynamics[END_REF] . The system is a cylindrical nanopore with radius R n = 2 nm and length L n = R n and connected to large half-spherical reservoirs with radius H n = 20 nm (Fig. 1 A B). The solid surfaces are made of a combination of fixed atoms and a solid flat wall that interacts with the liquid by generating a force on the atom in a direction perpendicular to the wall. Water is modeled only in an implicit way, as a medium that has the dielectric permittivity of water r = 80. The interaction between the solute and the water molecules is mimicked by applying a Langevin thermostat on the ions (temperature T = 300 K). The same damping time τ = 23 fs -controlling the amplitude of the viscous drag in the Langevin model -is imposed for both species, to ensure that they have the same diffusion coefficient, with a value D ≈ 1.5 • 10 -5 cm 2 /s, close to experimental ones for typical microions. The equations of motion are integrated using the velocity Verlet algorithm [START_REF] Martys | Velocity Verlet algorithm for dissipative-particle-dynamics-based models of suspensions[END_REF] , with a time-step equal to 2 fs. The atomic mass of the species is chosen to be m n+ = m n-= 39.1 g/mol, and we consider monovalent ions. A Weeks-Chandler-Andersen potential, keeping only the repulsive part of the Lennard-Jones potential, is used, and all particles (ions and solid atoms) share the same size σ = 3 Å. Long-range Coulombic interactions are calculated by means of a particle-particle particle-mesh (PPPM) solver.

Brownian particle simulations

We performed simulation of a random walk using a home-made c++ code. A number N of particle is initially introduced within a finite domain at random position. The domains are either a nanotube connected to half spherical reservoirs (Fig 1 A), an infinite cylinder