Karl Joulain

Analogy between thermal emission of nano objects and Hawking's radiation

We analyze in this work some analogies between thermal emission of nano objects and Hawking's radiation. We first focus on the famous expression of the black hole radiating temperature derived by Hawking in 1974 and consider the case of thermal emission of a small aperture made into a cavity (Ideal Blackbody). We show that an expression very similar to Hawking's temperature determines a temperature below which an aperture in a cavity cannot be considered as standard blackbody radiating like T 4 . Hawking's radiation therefore appear as a radiation at a typical wavelength which is of the size of the horizon radius. In a second part, we make the analogy between the emission of particle-anti particle pairs near the black hole horizon and the scattering and coupling of thermally populated evanescent waves by a nano objects. We show here again that a temperature similar to the Hawking temperature determines the regimes where the scattering occur or where it is negligible.

I. INTRODUCTION. HAWKING TEMPERATURE.

In 1974, Stephen Hawking [1,2] published his famous work about black hole radiating temperature. This famous temperature expression is the sole one in physics in which appears all 4 fundamentals constants (c, h, G and k b ). The origin of this radiation lies in the fact that the black hole makes vacuum fluctuations such as particle-antiparticle pairs radiate at infinity since they are separated at the black hole horizon : a particle is radiated to infinity whereas its antiparticle is pursuing its path into the black hole. The idea of the present work is two fold. First, we will remark that a similar formula for a critical temperature is obtained when one is studying thermal emission of an aperture. Second, an analogy is made between the particle-antiparticle pair separation at the horizon and the scattering of evanescent waves such as polariton by subwavelength structures. It is shown in this second part, that here again, this coupling to the far-field of evanescent waves thermally populated happens when the temperature is smaller than a critical temperature very similar to the Hawking one.

Hawking temperature formula can be heuristically derived from simple principles like the Heisenberg incertitude relation as well as simple relations between distance and velocity or wavelength and frequency. The idea is that in vacuum, there are quantum fluctuations that can be interpreted as permanent creation and annihilation of particle-antiparticle pairs. These pairs remain virtual as long as the Heisenberg incertitude relations are violated i.e the energy fluctuation ∆E due to the vacuum fluctuations and particle antiparticle pair formation multiplied by the time of their existence ∆t is smaller than the Planck constant (∆E∆t ≤ ). In the presence of a black hole horizon, the particle-antiparticle pairs close to the horizon can be separated. The typical time ∆t on which the pair are separated is the order of the time required to cross the whole horizon black hole that is ∆t ∼ 2πR s /c where R s = 2GM/c 2 is the Schwarzschild radius. Using the fact that the virtual pairs will be seen when the Heisenberg incertitude relations are valid, a particle that will be seen at infinite coming from the black hole will have their energy ∆E that will obey

∆E∆t ∼ /2 that is ∆E ∼ hc 8π 2 R s
T H = hc 3 16π 2 GM k b (1)
where the 4 fundamental constants appear. It is the temperature at which the black hole radiates.

II. THERMAL EMISSION BY AN APERTURE

In 2016, Joulain et al. [3] studied thermal emission by a subwavelength aperture. Indeed, a small aperture separating a cavity at thermal equilibrium from the outside is often seen as the archetype of a blackbody (Fig. 1). The idea is that such a structure is a perfect absorber since an incident radiation in the cavity will hardly be able to exit from. This remains true as long as the the aperture is large compare to the wavelength involved. When it is not the case, the aperture confines the radiation and diffracts it so that the transmission through the aperture is greatly reduced. This phenomenon is actually nothing else than the usual diffraction by a small aperture. In the case of a cavity which is at thermal equilibrium, the energy density in the cavity u(ω) is given by the so called Planck formula which is maximum at a frequency given by the Wien's law.

u(ω) = ω 2 π 2 c 3 (exp[ ω/k b T ] -1)
The angular frequency corresponding to the maximum of emission is given by

ω max = 2.8k b T
Therefore, one expects that as long as the aperture is much larger than the wavelength given by the Wien's law, it will behave and emit as a perfect blackbody. The condition reads in the case of an aperture of radius r 0

r 0 λ max = 2πc/ω max (2)
or if the temperature of the cavity is much larger than a critical temperature

T hc 2.8r 0 k b = T c (3)
Note that the right hand side of the preceding equation is very similar to the Hawking temperature except that some numerical factors are different and that the Schwarzschild radius is replaced by the aperture radius r 0 . Joulain et al. have solved the problem using fluctuationnal electrodynamics [START_REF] Rytov | Principle of Statistical Radiophysics 3, Elements of Radiation Fields[END_REF] when the aperture is larger that the wavelength (k 0 r 0 > 6 with k 0 = ω/c) or when k 0 r 0 1 (Bethe-Boukamp theory [START_REF] Bethe | [END_REF]6]). It is possible to choose at what precision an aperture can be assimilated to a blackbody. For example, if one plots the emissivity of such an aperture vs k 0 r 0 (See Fig. 2), one shows that this emissivity approaches 1 only for large k 0 r 0 . What we can emphasis here is that the aperture of a cavity at thermal equilibrium will behave like a blackbody really if its temperature is such that the photons emitted will have their wavelength smaller than the aperture. Similarly, it can be seen that a black hole emits radiation at a temperature which is such that the wavelength associated to this temperature is of the order of the Schwarzschild radius.

III. THERMAL EMISSION OF A SMALL STRUCTURE

In the last 20 years, many works have shown that thermal emission could be different whether it was considered in the near or in the far-field i.e at distances small or large compare to the thermal wavelength. For example, the density of electromagnetic energy of a single material and do not radiate in the far field. However, it is well known that these waves can be scattered by small objects, coupled to the far-field and can contribute to a strong increase of the radiation emission. This is the principle of scattering of surface waves by a grating [7] or by a microscopic tip [8]. The coupling to the far-field can be understood as followed. A surface wave has a parallel wave vector larger than k 0 so that the wave vector component perpendicular to the interface is pure imaginary in vacuum, evanescent and cannot propagate. The presence of grating, a tip or a discontinuity couple with the surface wave so that its parallel wave vector can be slightly modified in order that it goes below k 0 : the wave can now radiate into the far field. This phenomenon can also happen if the field is confined at a subwavelength distance. The analogy with black hole radiation remains here between the radiation of virtual particle-antiparticle pairs and between the scattering of the thermally populated surface waves by subwavelength objects.

In the case of an aperture filled with a material, Joulain et al. [3] gave an expression of the emitted flux which depends on a transmission function of the radiation through the aperture. When the transmission function is 1, one is in the situation of the Kirchhoff approximation so that the heat flux is given by

φ(ω) = φ 0 (ω) 2k 0 r 0 0 W (u/k 0 )uF (u)du × 1 0 κJ 0 (κu)dκ √ 1 -κ 2 (2 -|r s | 2 -|r p | 2 ) + ∞ 1 2κJ 0 (κu)dκ √ κ 2 -1 (r s ) + (2κ 2 -1) (r p ) e -2 √ κ 2 -1k 0 z (4)
where κ = K/k 0 , F (u) = (sin u -u cos u)/u 3 and

W (d) = 2 π arccos d 2r 0 -d 2r 0 1 -d 2r 0 2 .
W (u/k 0 ) is a function that decreases almost linearly from 1 to 0 when u is increasing from 0 to 2k 0 r 0 whereas F (u) decreases from 1 to 0 very fast as long as u is larger than 1. Under these conditions, if k 0 r 0 is much larger than 1, integration over u between 0 and 2k 0 r 0 can be replaced by an integration between 0 and ∞. As F (u) decreases fast with u when u is larger than 1, W can be approximated by 1 in the integral. Noting that

∞ 0 uF (u)J 0 (κu)du vanishes if κ > 1 and is equal to √ 1 -κ 2 if κ < 1 [9]
, one retrieves that there is no contribution of the evanescent waves to the emitted flux for large apertures and one retrieves in this case the classical expression of the emitted flux for a large surface

φ = φ clas = φ 0 (ω) 1 0 κdκ(2 -|r s | 2 -|r p | 2 ) = φ 0 (ω) dΩ cos θε(θ) (5)
where ε(θ) is the material emissivity, θ the angle with the normal to the surface and Ω the solid angle over which radiation angular integration is performed. From the preceding, discussion, one sees that the scattering of confined evanescent wave will occur when k 0 r 0 approaches 1. For thermal radiation, this corresponds to the condition where r 0 ∼ λ max i.e.

typical temperature of the order

T ∼ hc 2.8r 0 k b (6)
which is still very similar to the Hawking temperature where the Schwazschild radius has been replaced by the aperture size.

As in the case of vacuum aperture, this critical temperature corresponds to the one associated to the maximum radiation emission wavelength (thermal wavelength) at this temperature which is of the same order of the the object size.

This temperature also naturally appears when one is dealing with thermal emission of a perfect spherical absorber which size is much smaller than the thermal wavelength. Let us thus consider a sphere filled with a material which permittivity is . The scattering cross and absorption cross section can be calculated exactly with Mie Theory. Then, the total emission by the sphere is just the product of the emittance multiplied by the absorption cross section and integrated over all frequencies.

φ = ω 3 4π 2 c 2 (exp[ ω/k b T ] -1) σ abs dω (7)
In the case where the sphere size is much smaller than the thermal wavelength, the absorption cross-section is limited to the 1st term of the Mie Theory. The spherical object is in the dipolar approximation. Following Grigoriev et al. [START_REF] Grigoriev | [END_REF], the spherical object can exhibit a maximum cross section if its permittivity is suitably chosen. Note in that case that the absorption cross section is equal to the scattering cross section

σ max abs = σ max sca = 5π k 2 0 ( 8
)
Integrating over all the frequencies, the maximum flux emitted by a spherical object in the dipolar approximation is

φ max = 5πk 2 b T 2 24 = 1.18 × 10 -12 W K -2 T 2 (9)
which is apart to a numerical number the product of the quantum of thermal conductance with the temperature. Note that this flux does not depend on the sphere size (although the optimized dielectric function maximizing absorption cross-section σ abs depends on it). This

Note that this critical temperature is similar to the critical temperature T c that has been considered previously. Here again, when the temperature is such that the thermal wavelength associated is much smaller than the particle radius, the brightness temperature is larger than the actual temperature and is the geometric average of T and T c2 .

Finally, one can remember from the preceding considerations, that in analogy with black hole evaporation via the coupling of particle-antiparticle pairs, scattering of thermally populated evanescent waves occur at the surface of an object when their wavelength is of the order or smaller than the object size that is when its temperature is smaller than T c .

IV. CONCLUSION

We have shown in this work that a temperature very similar to the one introduced by

Stephen Hawking caracterizing the black hole evaporation can be introduced in the case of thermal emission of a small object. When the thermal wavelength associated to thermal emission is smaller than the size of the object, it is known that this thermal emission greatly differs from the one given by usual radiometry i.e geometrical optics. This allows to introduce a critical temperature very similar to the one introduced by Stephen Hawking for black hole radiation. Moreover, when the preceding conditions are reached i.e when the temperature is smaller than the critical temperature, an analogy can be done between the separation of particle-antiparticle pairs due to vacuum fluctuations by a black hole horizon and the scattering of evanescent surface waves due to thermal fluctuations by nano objects.

FIG. 1 :

 1 FIG. 1: Cavity at thermal equilibrium at temperature T separated from the outside by an aperture of radius r 0 .

FIG. 2 :

 2 FIG. 2: Emissivity of the cavity versus k 0 r 0 calculated by means of fluctuation electrodynamics in the Bethe-Bouwkamp regime (k 0 r 0 1, dashed-green curve) and in the Kirchhoff regime (k 0 r 0 > 6, blue curve).

σT 4 b ( 10 )

 410 maximum flux can be compared to the classical emitted flux that would be given in the case of a perfect spherical blackbody of radius r 0 with brightness temperature T b φ clas = 4πr 2 0 where σ = π 2 k 4 b /60 3 c 2 is the Stefan-Boltzmann constant. The brightness temperature therefore reads T b = √ T c2 T where


Julien Delplanque
email: julien.delplanque@inria.fr

Stéphane Ducasse
email: stephane.ducasse@inria.fr

Andrew P Black
email: apblack@pdx.edu

Guillermo Polito
email: guillermo.polito@inria.fr

Guillermo Black

Polito

Rotten Green Tests A First Analysis

Keywords: rotten green tests, unit tests, testing SetTest >> testAddTheSameElementTwiceResultOneOccurrence

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Agile methodologies such as Extreme Programming [START_REF] Beck | Manifesto for Agile Software Development[END_REF][START_REF] Beck | Test Driven Development: By Example[END_REF][START_REF] Beck | Extreme Programming Explained: Embrace Change (2Nd Edition)[END_REF] promote Unit Testing [START_REF] Meszaros | The Test Automation Manifesto[END_REF] as a key tenant of the software development process. Executing a test suite after each change to the software helps to ensure that new functionality works, and that the old functionality remains working, that is, it helps avoid software regressions [START_REF] Astels | Test-Driven Development -A Practical Guide[END_REF].

Tests are based on the execution of assertions that check that the system under test satisfies some property, for example, that a method returns a certain value, or that certain data is written to a stream. Developers value "green tests", i.e., tests that are passing, because they provide assurance that the software is working as expected.

RMOD Inria Technical Report, Jun -2018, 2018.

Our concern in this work is with tests that were intended by their designer to execute some assertions, but do not actually do so -we call them rotten green tests. Such tests are insidious because they pass, and they contain assertions; they therefore give the impression that some useful property is being validated. In fact, a rotten green test guarantees nothing: it is worse than having no test at all! Our approach is based on a combination of static analysis and dynamic monitoring of method execution. We identify whether or not a test is rotten, even in the presence of helper methods and trait compositions. A limitation of our current implementation is that a test with multiple assertions will be considered as rotten only when none of the assertions is executed. Future work will ensure that all the assertions are executed.

The contributions of the paper are:

• the recognition of rotten "green" unit tests that contain assertions that are not executed, and which therefore give the developers false confidence (Section 2); • a simple combination of static and dynamic analyses that identifies such rotten green tests, (Section 3); and • a report on our experience applying this approach to several large systems (Section 4).

The Problem of Rotten Green Tests

Before defining rotten green tests, we first describe the basics of unit testing, and then briefly explain "Smoke Tests", to help distinguish them from the topic of this article.

Unit tests

Unit tests are commonly composed of a test fixture (which sets up the system to be tested), one or more stimuli (which exercise the component under test), and one or more assertions that verify some expected property [START_REF] Black | Pharo by Example[END_REF][START_REF] Meszaros | XUnit Test Patterns -Refactoring Test Code[END_REF]. The following trivial example shows an SUnit test that checks that a set should not contain the same object twice.

In the example, the fixture is the code that declares and initializes s to contain 1; here the fixture is inline, but it can also be factored-out into a setUp method. The stimulus is the second addition of 1 to s; the assertions then verify the property that s contains 1 just once.

Provided that all of the assertions are true, this test will pass; we say that it is "green". If a false assertion is executed, for example, if the set does not detect the duplicate and its size is 2, the test will fail: it will be "yellow". If an error occurs during the running of the test, for example, if Set new signals an exception, then the test will be "red".

Smoke Tests

It is common practice to use unit testing frameworks to execute so-called "smoke tests" whose purpose is to check that the feature under test can be run without emitting "blue smoke" -that is, the test ran without raising an unexpected exception [START_REF] Waletzky | Smoke Tests vs[END_REF]. Here is an example of such a test.

SetTest >> testSetAddSmokeTest

| s | "Fixture" s := Set new. s add: 1.

"Stimulus" s add: 1

In its simplest form, a smoke test may contain no assertions at all, as in this illustration. This is the way that we use the term smoke test in the remainder of this article.

Smoke tests are useful because, if they are green, they provide a fast but cursory check that the feature concerned can be considered for further testing. Conversely, if a smoke test is red, there is a serious issue that should be addressed rapidly. Smoke tests are not the concern of this article; nothing that follows should be construed as advocating either for or against the use of smoke tests. Nevertheless, we do need to distinguish a smoke test that by design contains no assertions, from a rotten green test, which by accident executes no assertions.

Rotten Green Tests

Consider an empty test -a test method that contains no code at all: no fixture, no stimulus, and no assertion. If it is treated as a passing test, it will increase the number of green tests A much more insidious problem is caused by a test that does contain a valid fixture, stimulus and assertion, but which nevertheless does not execute any assertions. How can this happen? Let's look at a real example, taken from Pharo issue 7478, and shown in Figure 1.

At first glance, this looks like a fine test of the collection printOn:delimiter: method. The first three lines of the test method create a fixture -in this case they set up aStream. Then comes the stimulus: the collection self nonEmpty is sent the printOn: aStream delimiter: ',' message, causing it (we hope) to write its elements to aStream separated by commas. The remainder of the test is intended to make assertions. The code looks as though it is parsing the contents of aStream and asserting that the elements it finds written there are the same as those in the original collection.

This test is green, so we know that everything is OK, right? Wrong! The programmer who wrote this test misunderstood the way that streams work. The string result can never be modified by writing to aStream; indeed, the very name result is misleading, because result is initialized to the empty string and never changes. Consequently, result find-BetweenSubstrings: will answer an empty collection, and the withIndexDo: block -which contains the only assertion -will never be executed. We could put any assertion into this block, and the test would still run green.

We believe that such a test is worse than no test at all. First, the rotten test assures us of no property of the method under test. Second, it wastes time -remember that we want our tests to run quickly. Third, in the case of "no test at all", test coverage statistics might reveal that the printOn:delimiter: method is not being exercised, and the developers would at least be aware that their testing is inadequate. Instead, with a rotten test in place coverage will tell us that the printOn:delimiter: method is being executed -and we will be unaware that it is not actually being tested.

Rotten green tests give developers a false sense of security: coverage may be good, and the tests may be green, but in fact no properties are being checked (other than the property that the code runs without signaling an error). Moreover, when developers look at such a test superficially, they will probably not spot that there is a problem.

The example in Figure 1 was first reported as a bug in February 2013. However, because the test was green, the bug was easy to overlook, and the bug report was closed without any action being taken. As of the end of 2017, this rotten green test was still present in Pharo.

Identifying Rotten Green Tests

Once we accept that rotten green tests are bad, it is natural to ask how we can detect them. One might think that all that is necessary is to detect tests that makes no assertions, but this won't let us distinguish smoke tests from rotten green tests. Can't we identify smoke tests as those tests that contain no assertions? No, because many tests make assertions indirectly through the use of helper methods (as explained in Section 3.1); tests that use assertions in helper methods are not smoke tests, even though they contain no assertions directly.

To clarify the discussion, we will use the following terms.

• an assertion primitive is a method of the unit-testing framework that performs the actual check. In Pharo, assert: and assert:description: are the only assertion primitives. • a test method is a method identified as such by the unittesting framework. In Pharo, test methods are zeroargument methods defined in a subclass of TestCase whose names start with "test". • a helper method is a method that makes an assertion directly (by invoking an assertion primitive) or indirectly (by invoking another helper method), but that is not a test method. In Pharo, SUnit provides helper methods like assert:shouldRaise:; in addition, developers frequently write their own application-specific test helper methods, as we will discuss in Section 3.1.

In Figure 2, testABC is a test method, and helper and sec-ondHelper are helper methods, because helper invokes sec-ondHelper, and secondHelper invokes an assertion primitive.

Helper methods

It is common practice for developers to factor-out assertions into helper methods. This affects our analysis in two ways: we need to know which messages invoke helper methods when we identify smoke tests, and we must take into account the possibility that helper methods, as well as tests, might be rotten.

What is a helper method? In our context, it is a method that makes an assertion, either using another helper method, or by using one of the assertion primitives. As a simple example, a test suite for approximate numerical methods, might define and use the helper method assert:isRoughly:within: Like a rotten test, a helper method might fail to make an assertion in some or all situations; in other words, helper methods might be rotten too. Any approach to detecting rotten green tests should also detect rotten helper methods. Moreover, since the action of a helper method may depend on the context (e.g., the test fixture, and the arguments to the helper method), a helper might work fine in one context, but be rotten in another. So, when detecting the rotten helper, we need to record the specific test that exposes the rot.

In Pharo, traits enable developers to reuse tests across several sub-hierarchies; the Collection tests are a good example [START_REF] Ducasse | Traits: A Mechanism for fine-grained Reuse[END_REF][START_REF] Ducasse | Reusing and Composing Tests with Traits[END_REF]. The situation with traits is the same as with helper methods: a trait method may show up as being rotten only for a certain trait use. This can happen because trait composition and inheritance can change the test fixture.

Classifying Tests

There are three situations that a rotten test analysis should identify.

Good tests. A test passes, contains some assertions (either directly, or indirectly though helper methods), and some assertions are executed: the test is good. Rotten tests. A test passes, and contains assertions (either directly, or indirectly), but no assertion is executed: the test, or the helper method, is rotten. Smoke tests. A test contains no assertions (either directly, or indirectly); it is a smoke test.

Note that distinguishing between good and rotten requires some dynamic analysis, because we need to ascertain whether an assertion is executed. In contrast, distinguishing between rotten and smoke requires some static analysis, because we need to ascertain whether the test contains assertions.

Combining Static and Dynamic Analyses

Our analysis performs the following steps:

Step 1: identification of assertion primitives. We build the set of assertion primitives by identifying all the methods of the unit test framework that make assertions directly. This set depends only on the test framework.

Step 2: identification of helper methods. To build the set of helper methods, we start with a set S containing all the methods in the test class that are not test methods but which self-send one of the assertion primitives. ("All the methods" includes methods introduced through inheritance or trait use.) We identify such methods using a simple static analysis of self-and super-sends. We then add to S any method in the test class that is not a test case but which self-or super-sends one of the methods in S. We repeat this step until no new methods are added to S; the set S now contains all the helper methods.

Step 3: test execution. We execute each test method (including inherited test methods) one at at time, while monitoring the following pieces of information: a. the outcome of the test (pass, fail, or error), b. whether or not one of the helper methods is executed, and c. whether or not one of the assertion primitives is executed. Because we are looking for rotten green tests, we consider only passing tests.

Step 4: classification. We combine the test execution information from Step 3 with the static information collected in Steps 1 and 2. Table 1 shows how this information is used to classify test methods. We discuss the classification in Section 3.4.

Step 5: report generation. DrTest's final report has to take into account the way that methods are reused in the test hierarchy. A test method may be defined in a superclass and executed in a subclass. A test method may also be defined in a trait [START_REF] Ducasse | Traits: A Mechanism for fine-grained Reuse[END_REF], and reused by several classes. In general, the test fixture, and thus the meaning of the test, will be different in each place in which it is used, so we must report the class of the test as well as the method. Helper methods are designed to be used by many test methods; if a helper method is rotten, the test invoking it should be reported, so that we can understand the scenario in which the helper fails to make an assertion.

Our current dynamic analysis has a coarse granularity: we monitor the execution of primitive assertions for the complete run of the test method under analysis. This means that we cannot attribute an executed assertion to a specific call site. So, if one primitive assertion in the test method was executed, the test will not be classified as rotten, even though there might be another primitive assertion that was not executed. We discuss this shortcoming further in Section 6.

Classification Discussion

We now discuss the reasoning behind the classification in Table 1. In the first four rows, the test contains (statically) sends of both assertion primitives and helper methods. If the test is good, we will see both helper and assertion executed, as in row 1. If no helper, or no assertion primitive, or neither, is executed, then the test is rotten.

In rows 5-8, there is no send to a helper method in the test, but there is a send of an assertion primitive. If an assertion is executed (rows 5 and 6), the test is good; if no assertion is executed (rows 7 and 8), we can say that the test is rotten. On rows 5 and 7, a helper method is executed even though there is no send of a helper message in the test! How can this happen? One possibility is that the test constructs the selector of the helper method dynamically, and then uses perform: to send it. Another possibility is that the test invokes a helper method with a message send that is not a self-or super-send, either because the message is actually sent to an object outside its class hierarchy, or because the receiver of the message is not the pseudo-variable self. We flag these entries with "dynamic helper invocation"; test like this are unlikely to occur in practice, but if they do, they are worth a second look.

Rows 9-12 are similar to rows 1-4, except that in rows 9-12 all of the assertions are made using helper methods. If both helper and assertion are executed (row 9), the test is good; if either is not executed, something is rotten. In row 11, the helper is executed, but the assertion is not, so it is the helper method that is rotten. In rows 10 and 12, no helper is executed, even though the test relies on helper methods, so we know that the test is rotten. This reasoning applies even in row 10 where, somehow, an assertion is executed even though no assertion message is sent!

In rows 13-16, no helpers and no assertions are used in the test. Such tests look like smoke tests -in which case we would expect that no helpers and no assertions are executed, 2. The dynamic analysis will report that neither a helper method nor an assertion primitive is executed. The static analysis will identify secondHelper as a helper method (because it makes a primitive assertion), and helper as a helper method (because it uses secondHelper). Hence, we look for ✗✗✓✗ in Table 1, which we find on row 12. As expected, this tells us that we have found a rotten test.

Results

Pharo is an open-source language with a growing community and a large number of mature projects. We have run DrTest on eight Pharo subsystems; we found that all subsystems except Iceberg have at least 1 rotten test. The results are shown in Table 2 Calypso is the new system browser for Pharo. It allows developers to view and edit packages, classes, and methods. DrTest found four rotten tests; two of them contained conditionally-executed sends of TestAsserter » assert:description: with false as the first argument, as shown in Figure 3. These two tests are false positives: the developer intends that the assertion will not be executed if the code under test behaves correctly, that is, if it signals an exception. This test would be clearer if it simply asserted that ClyCompositeScope on: #() should raise an error.

The two other rotten tests in Calypso contain helper methods that are never executed. These tests are tagged with the <expectedFailure> pragma which indicates that they are expected to fail by SUnit. Again, those are false positives.

Andrew ▶I don't see why, but then, I can't find the tests.◀ Collections are the basic data structures of Pharo [START_REF] Black | Pharo by Example[END_REF]. The collections tests contain seven rotten tests. Two of them are tests that contain a guard clause checking which of the collection classes is being tested. This guard clause is present because the test is implemented in a superclass common to multiple test cases, but the test should be run only on certain subclasses. One rotten test contains a guard clause that deliberately skips the test when the Dictionary subclass under test does not support nil keys.

The remaining four rotten tests are true positives that use the Stream API incorrectly. These methods test the serialization of a collection on a stream, using the following process: (1) create a stream, (2) write a collection on the stream and ( 3) compare what was written on the stream to the original collection. As discussed in Section 2.3, the problem here is that the comparison should be with the contents of the stream, aStream contents, and not with the variable result, which remains bound to the empty string. Because of this, the assertion inside the loop is never executed. Figure 4 illustrates another such rotten test.

Glamour provides a high-level API for creating user interfaces. It has three rotten tests. After an investigation, we realised that two of them are expected failures. Andrew ▶So what? If the tests fail, then they can't be rotten, so I assume that they pass. Are they rotten?◀ The last one is because a block containing an assertion is never executed Andrew ▶Is this a true positive? Lets say so◀.

Iceberg is a tool to manage Git project within Pharo. It was found to contain no rotten tests.

OpalCompiler is the default compiler used in Pharo images. It has 15 rotten tests; 10 of them are due to a guard clause that was incorrect. As shown in Figure 5, the comment above the guard clause specifies that the test should not be run unless the method's source code is available. But the implementation of the condition in the guard clause ensures that the test is skipped if the source code is available.

The other five rotten tests of the OpalCompiler were intended to identify a bug related to dynamic bytecode rewriting of boolean expressions. In Smalltalk, it is normally impossible to make a non-Boolean object act as a Boolean in a condition, even if the object implements the Boolean interface. The reason is that common Boolean message like #ifTrue: and #and: are compiled to optimized bytecode that raises an exception when they are sent to non-Booleans. However, Pharo dynamically catches this exception and rewrites the bytecode to a de-optimized version. This allows one to use a non-Boolean receiver as a Boolean (provided that it implements the methods of Boolean). The five rotten OpalCompiler tests are intended to validate this feature. A bug in the bytecode rewrite process introduced an early return into the test methods. The return appears before the assertion can be executed. Thus, it leads to five tests passing but not executing any assertion. Figure 6 shows the source code of one of these five methods, both as it appears, and as it is dynamically rewritten. These rotten tests teach us that it can be extremely difficult to see whether or not an assertion is executed.

Pillar is a markup syntax and associated tools for writing and generating documentation, books, and slides. This project contains 1 rotten test which is tagged as an expected failure. Andrew ▶again, what has expected failure to do with things?◀ System packages are also provided with test suites. These test suites contain one rotten test, which is tagged as an expected failure. Andrew ▶again, what has expected failure to do with things?◀ Zinc is a library for interacting with the network. It contains three rotten tests. One of them is a test that sends #skip, a method of TestAsserter whose specification says "Don't run this test, and don't mark it as failure". Thus, we can be confident that the assertions in this test are skipped on purpose.

Another Zinc test is guarded by a check that a certain class exists in the system; if it does not, the test returns. By default, the class does not exist, so the test returns and no assertion is executed.

The last rotten test of Zinc performs a network HTTP operation. If the result is a redirection, the test returns immediately. The URL requested by the test does lead to a redirection, so the test returns and no assertion is executed.

Implementation

The implementation of DrTest has two points of interest: the detection of helper methods and the detection of method execution. Helper methods are detected through static analysis: we traverse the class hierarchy starting from the subclass of TestCase that interests us, and classify as helper methods all non-test methods that make a direct or indirect send of an assertion primitive. For the dynamic analysis, the detection of method execution is done by instrumenting the assertion primitives and helper methods. We replace those methods by method spies that record that they were called, and then forward execution to the original method.

The Algorithm in Outline

To identify the rotten tests of a test class, DrTest starts by determining which methods are helpers (Section 5.2) and marking them with spy objects containing a flag called (Section 5.3). Then, DrTest executes each test case by:

1. reseting all spies so that called = false, 2. running the test, which will mark executed method by setting called = true, and 3. classifying the outcome of the test using the marked spies according to Table 1. The code that classifies a test outcome is as follows.

• Table Rows

Detecting Helper Methods

DrTest detects helper methods through static analysis. The analysis first collects all the test methods Andrew ▶do we mean the non-test methods?◀ in the test class under analysis, and all of its superclasses up to TestAsserter (which is the superclass of TestCase). All of these methods may be relevant to the execution of the test class under analysis. We then use an AST visitor to visit the entire AST of each of these methods; every time the visitation encounters a self-send, it notes the method corresponding to the self-send. This gives us an enumeration of the relation "executes directly via a self-send". We then compute the transitive closure of this relation, and note the methods that directly or indirectly self-send an assertion primitive. The noted methods are the helper methods.

Detecting Method Execution

The method spy mechanism uses the method wrapper virtual machine (VM) hook described by Martinez Peck et al. [START_REF] Peck | Ghost: A Uniform and General-Purpose Proxy Implementation[END_REF], which works as follows. When a message is sent to an object, the VM looks in the method dictionary of the receiver's class for an entry keyed by the method selector. In most cases, this entry will be a compiled method, that is, an instance of the class CompiledMethod, in which case the VM executes the method with the arguments from the message. However, when the object found in the method dictionary is not a compiled method, the VM instead sends the message run:with:in: to the found object, giving it a chance to execute some behaviour. This hook is commonly used by profilers and other dynamic analysis tools.

In our case, the object in the method dictionary is an instance of MethodTracer, a class that implements a run:with:in: method that contains a called flag, which is initialized to false, and a reference to the original method. When a Method-Tracer is sent a run:with:in: it first "marks" itself by setting the called flag to true, and then forwards execution to the original compiled method. Here is the relevant code from

MethodTracer:

MethodTracer >> run: aSelector with: someArguments in: aReceiver self mark. ^aReceiver withArgs: someArguments executeMethod: method Distinct instances of MethodTracer are used to spy on the assertion primitives and on the helper methods. This allows us to keep a separate record of whether or not assertion primitives and helper methods have been executed.

Discussion and Future Work

As currently implemented, DrTest has two significant limitations.

Tests containing multiple assertions

When a test method uses multiple assertion primitives or helper methods, our approach does not currently distinguish between the case in which just one assertion primitive or helper methods is executed, and the case where they are all executed. In both cases, we report that the test is good.

For example, the test in Figure 1 is rotten, and is detected as such. However, if the developer had added self assert: (self nonEmpty isNotEmpty) description: 'test fixture broken' at the top of the test, it would not be detected as rotten, even though this assertion does nothing to solve the problem. Indeed, because it is based on method-level monitoring, the current implementation cannot distinguish between assertions made at different call sites. This can result in false-negatives, that is, not all rotten tests will be detected. It will not cause false-positives: a test that is identified as rotten will always be rotten.

We cannot yet evaluate the number of false negatives due to multiple assertions, because doing so would require that we enhance our solution to eliminate them.

For the future, we can see several ways of dealing test that contain with multiple assertions. One possibility is to create, from each test method that contains multiple assertions, multiple tests, each with a single assertion, and then to apply our process to these new methods. Another possibility is to record, not just when an assertion is executed, but which send of the assertion message is being executed. These enhancements are the subject of ongoing work.

Location of Helper Methods

Another limitation of our approach is that our static analysis discovers only those helper methods that are in the hierarchy of the analysed test case. If the helper methods are located in a utility class, then they will not be detected. We can estimate how frequently this problem arises because some instances will manifest themselves as a test method executing an assertion even though it does not contain any assertions or helper methods -the cases classified as involving a dynamic invocation in Table 1. (Tests that contain self-sends that execute helper methods, and sends of helper messages to a utility object, will not be detected.) Our preliminary results indicate that these cases do not occur frequently.

Related Work

Software testing is an active area of research; many researchers have looked at improving the quality of tests, but we are not aware of any prior work that identifies rotten green tests.

Mutation testing is one of the earliest approaches used to improve test quality and robustness [START_REF] Demillo | Hints on Test Data Selection: Help for the Practicing Programmer[END_REF]. Mutation testing generates "mutant" programs from the program under test, and then selects tests that differentiate the mutants from the original program. Tests that detect few or no mutants are candidates for removal.

Several researchers have used mutation testing to improve branch coverage [START_REF] Lingampally | A Multipurpose Code Coverage Tool for Java[END_REF]. Tillmann and Schulte [START_REF] Tillmann | Parameterized unit tests[END_REF] use symbolic execution to find inputs for parameterized unit tests that achieve high code coverage. They turn existing unit tests into parameterized unit tests and generate entirely new parameterized unit tests that describe the behaviour of an existing implementation.

Baudry et al. [START_REF] Benoit Baudry | Automatic Test Case Optimization: A Bacteriologic Algorithm[END_REF] present a bacteriological approach to mutation testing.

Other approaches focused on other attributes of test quality. Baudry et al. [START_REF] Benoit Baudry | Improving test suites for efficient fault localization[END_REF] addresses improving the value of tests for diagnosis. They propose a new attribute called the Dynamic Basic Block, to assist in the task of locating faults in a program. For fault localization, the usual assumption is that test cases satisfying a chosen test-adequacy criterion are sufficient to perform diagnosis. This assumption is verified neither by specific experiments nor by intuitive considerations.

Often, unit test frameworks present failed unit tests in an arbitrary order, but developers want to focus on the most specific ones first. Gaelli et al. [START_REF] Gaelli | Ordering Broken Unit Tests for Focused Debugging[END_REF] propose a partial order of unit tests corresponding to the containment hierarchy of their sets of covered method signatures. When several unit tests fail in the same hierarchy, the tool can guide the developer to the test involving the smallest number of methods.

Other work focuses on the selection of the tests to be run. For example, when a change is made to some software, it is desirable to rerun those tests that are most likely to be invalidated by the change. Beszedes et al. [START_REF] Beszedes | Code Coverage-based Regression Test Selection and Prioritization in WebKit[END_REF] propose to use code coverage for test selection testing to maximize the test surface. Blondeau et al. [START_REF] Blondeau | Test Case Selection in Industry: An Analysis of Issues Related to Static Approaches[END_REF] analyse the problem of test selection surfaces in an industrial context. Verhaeghe et al. [START_REF] Verhaeghe | Usage of Tests in an Open-Source Community[END_REF] examined the practices of developers when writing code, and the way they execute (or choose not to execute) tests. SmartTest [START_REF] Demeyer | Evaluating the Efficiency of Continuous Testing during Test-Driven Development[END_REF][START_REF] Verhaeghe | Usage of Tests in an Open-Source Community[END_REF] is a tool that automatically selects the tests to be run.

Deursen et al. [START_REF] Van Deursen | Refactoring Test Code[END_REF] present a list of "bad test smells" and their associated cures. They do not mention rotten green tests as a smell.

The term "smoke test" is used in various ways by different authors. Waletzky [START_REF] Waletzky | Smoke Tests vs[END_REF] says that the terms Smoke Test and Build Verification Test are sometimes used interchangeably, but prefers to treat Smoke Tests as the subset of Build Verification Tests that are extremely fast to run, and are the prelude to more thorough testing. Other authors use the term "smoke test" to include tests that make assertions. For example, Memon and Xie [START_REF] Atif | Empirical Evaluation of the Fault-Detection Effectiveness of Smoke Regression Test Cases for GUI-Based Software[END_REF] discuss generating thousands of tests that contain sequences of simulated GUI events, and using various test oracles to check that the state of the GUI is as expected -which they still call smoke tests.

The implementation of DrTest uses a form of method wrapper to count executions of helper methods and assertions. Brant et al. [START_REF] Brant | Wrappers to the Rescue[END_REF] discuss many techniques for implementing wrappers. Our approach relies on a Pharo VM hook that was not available to Brant et al., who instead subclassed compiled method. Martinez Peck et al. [START_REF] Peck | Ghost: A Uniform and General-Purpose Proxy Implementation[END_REF] describe the VM hook used by DrTest; they use it to implement Ghost proxies.

Conclusion

We have identified the existence of Rotten Green Tests, that is, tests that pass and contain assertions, but whose assertions are not executed. Such tests are worse than no tests at all, because they give developers false confidence in the system under tests. We have described an algorithm that identifies Rotten Green Tests, based on a combination of static and dynamic analysis, and a tool, DrTest, that implements the proposed approach. DrTest distinguishes rotten green tests from smoke tests (which also execute no assertions, but do so by design). We report on the tests found in eight large open-source projects.

DrTest is publicly available on github and can be loaded from github://juliendelplanque/RottenTestsFinder/src.

Figure 3 .

 3 Figure 3. Rotten test in Calypso. The assertion is not reached on purpose, this is a false positive.

  TPrintTest >> testPrintElementsOn | aStream result allElementsAsString tmp | result:=''. aStream:= ReadWriteStream on: result. tmp:= OrderedCollection new. self nonEmpty do: [:each | tmp add: each asString]. self nonEmpty printElementsOn: aStream . allElementsAsString:=(result findBetweenSubstrings: ' ' ). 1 to: allElementsAsString size do: [:i | self assert: (tmp occurrencesOf:(allElementsAsString at:i)) = (allElementsAsString occurrencesOf:(allElementsAsString at:i)). ].

Figure 4 .

 4 Figure 4. A Rotten test of Pharo Collections. The stream API is used incorrectly.

Figure 5 .Figure 6 .

 56 Figure 5. Rotten test in OpalCompiler. The guard clause implements the negation of what is described in the comment.

  Figure 2. A rotten test that uses a helper method.

	RottenTest >> testABC
	"Test method"
	false ifTrue: [self helper]
	RottenTest >> helper
	"Indirect helper"
	self secondHelper
	RottenTest >> secondHelper
	"Direct helper"
	self assert: x

NumericalTests >> assert: actual isRoughly: desired within: tolerance self assert: (actual -desired) abs <= tolerance description: [ 'actual result ' , actual , ' is not even roughly equal to ' , desired ] For a more realistic example, we turn to the Pillar editing platform [11]. Pillar's test suites use helper methods such as assertWriting:includesText:; which is defined in the superclass of Pillar's test classes, PRDocumentWriterTest. This method factors out the creation of a PRAnchor and its verification. Here is its definition and a sample of its use: PRDocumentWriterTest >> assertWriting: anItem includesText: aString | result | result := self write: anItem. self assert: result includesSubstring: aString PRHTMLWriterTest >> testAnchor | item | item := PRAnchor new name: 'foo'. self assertWriting: item includesText: 'id="foo"'

Table 1 .

 1 Classifying tests using the results of our analyses. In all cases, the test passes. "Executed" means that the (helper or assertion) method was executed. "Contains" means that the test method contains (statically) a message send to a helper or to an assertion.

		Dynamic Analysis	Static Analysis
	Row	Helper	Assertion Test contains Test contains
	№	Executed Executed	helper	assertion	Classification
	1	✓	✓	✓	✓	✓Good test
	2	✗	✓	✓	✓	✗ Rotten test
	3	✓	✗	✓	✓	✗ Rotten test & rotten helper
	4	✗	✗	✓	✓	✗ Rotten test
	5	✓	✓	✗	✓	✓Good test (dynamic helper invocation)
	6	✗	✓	✗	✓	✓Good test
	7	✓	✗	✗	✓	✗ Rotten test & rotten helper (dynamic helper invocation)
	8	✗	✗	✗	✓	✗ Rotten test
	9	✓	✓	✓	✗	✓Good test
	10	✗	✓	✓	✗	✗ Rotten test (dynamic assertion invocation)
	11	✓	✗	✓	✗	✗ Rotten helper
	12	✗	✗	✓	✗	✗ Rotten test
	13	✓	✓	✗	✗	✓Good test (dynamic assertion & helper)
	14	✗	✓	✗	✗	✓Good test (dynamic assertion invocation)
	15	✓	✗	✗	✗	✓Good test (dynamic helper invocation)
	16	✗	✗	✗	✗	✓Smoke test

as occurs on row 16. If a helper or assertion is executed, it is likely to be because of a dynamic invocation.

Let us see how this classification works on the example in Figure

Table 2 .

 2 . Rotten tests in Pharo subsytems.

	Subsystem	Packages Classes Test classes Tests Rotten tests
	Calypso	58	705	128 2671	4
	Collections	16	224	59 5858	7
	Glamour	19	463	65	449	3
	Iceberg	16	565	44	555	0
	Opal Compiler	7	227	49	854	15
	Pillar	33	358	112 3188	1
	System	48	330	44	552	1
	Zinc	9	184	43	412	3

ClyCompositeScopeTests>>#testEmptySubscopesAreForbidden

[ClyCompositeScope on: #(). self assert: false description: 'empty subscopes should be forbidden'] ifError: [].

  1-4.

	(containsAssertionPrimitive not and: [
	containsHelper and: [
	assertionPrimitiveExecuted not and: [
	helperExecuted ] ] ])
	ifTrue: [ self addRottenHelper: compiledMethod ].
	• Table Rows 10 and 12.
	(containsAssertionPrimitive not and: [
	containsHelper and: [
	helperExecuted not ] ])
	ifTrue: [ self addRottenTest: compiledMethod ].
	(containsAssertionPrimitive and: [
	containsHelper and: [
	(assertionPrimitiveExecuted and: [
	helperExecuted ]) not ] ])
	ifTrue: [ self addRottenTest: compiledMethod ].
	• Table Rows 7-8.

(containsAssertionPrimitive and: [ (containsHelper not) and: [ assertionPrimitiveExecuted not ] ]) ifTrue: [ self addRottenTest: compiledMethod ]. • Table Row 11.

Acknowledgements

This work was supported by Ministry of Higher Education and Research, Nord-Pas de Calais Regional Council, CPER Nord-Pas de Calais/FEDER DATA Advanced data science and technologies 2015-2020. Andrew Black was able to spend time at Inria because of a sabbatical granted by Portland State University.


R Lucken

V Croes

T Lafleur

J.-L Raimbault

A Bourdon

P Chabert

Edge-to-center plasma density ratios in two-dimensional plasma discharges

Keywords: Pre-sheath drop, edge-to-center plasma density ratio, low-temperature plasma discharge

Edge-to-center plasma density ratios -so-called h factors -are important parameters for global models of plasma discharges as they are used to calculate the plasma losses at the reactor walls. There are well-established theories for h factors in the one-dimensional case. The purpose of this paper is to establish h factors in twodimensional (2D) systems, with guidance from a 2D particle-in-cell (PIC) simulation. We derive analytical solutions of a 2D fluid theory that includes the effect of ion inertia, but assumes a constant (independent of space) ion collision frequency (using an average ion velocity) across the discharge. Predicted h factors from this 2D fluid theory have the same order of magnitude and the same trends as the PIC simulations when the average ion velocity used in the collision frequency is set equal to the ion thermal velocity. The best agreement is obtained when the average ion velocity varies with pressure (but remains independent of space), going from half the Bohm velocity at low pressure, to the thermal velocity at high pressure. The analysis also shows that a simple correction of the widely-used 1D heuristic formula may be proposed to accurately incorporate 2D effects.

Introduction

Global models are handy tools to assess the scaling laws of a plasma discharge with respect to parameters such as the background pressure, the input power, and the geometry in numerous applications, for instance to predict the performances of plasma thrusters [START_REF] Goebel | Fundamentals of Electric Propulsion: Ion and Hall Thrusters[END_REF][START_REF] Grondein | Global model of an iodine gridded plasma thruster[END_REF][START_REF] Chabert | Global model of a gridded-ion thruster powered by a radiofrequency inductive coil[END_REF]. They are zero-dimensional (0D) models that yield an approximate solution of particle balance and power balance equations [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF][START_REF] Chabert | Radio-Frequency Plasmas[END_REF]. It was shown in the past that these models are in good agreement with experiments [START_REF] Perret | [END_REF], as well as with more detailed fluid or particle-in-cell (PIC) simulations [START_REF] Lafleur | Edge-to-center density ratios in low-temperature plasmas[END_REF].

In low temperature plasma discharges, the region near the wall is characterized by a sheath where quasi-neutrality is violated. According to the Bohm sheath criterion, ions must leave the bulk plasma and enter the sheath region at a velocity higher or equal to the Bohm velocity u B = kTe

m i 1/2
, T e being the electron temperature, m i the ion mass, and k the Boltzmann constant. Since ionization is usually negligible in the sheath (at least in the low-pressure regimes considered in this paper), the ion flux reaching the wall is:

Γ i = n s u B ( 1
)
where n s is the plasma density at the sheath edge. For global models however, the flux at the wall should be expressed as a function of the space-averaged density, which is not very different from the density at the center, noted n 0 . It is then useful to introduce an h factor, defined as

h = n s /n 0 (2)
In the literature, this factor has been labeled h L for 1D Cartesian geometries or h R for 1D cylindrical geometries.

Several studies provided a good estimate of these factors over large pressure ranges, using various simplifications of the momentum conservation equations in one dimension (1D) [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF][START_REF] Chabert | Radio-Frequency Plasmas[END_REF]. The transport in the higher pressure regime was first investigated by Schottky [START_REF] Schottky | Diffusionstheorie der positiven Säule[END_REF]. In this limit, the ion-neutral collision frequency only depends on the ion thermal velocity

v i = 8kT i πm i 1/2
, where T i is the temperature of the ions. The lower pressure (collisionless) regime was studied by Tonks and Langmuir in the 1920's [START_REF] Tonks | A General Theory of the Plasma of an Arc[END_REF]. The intermediate pressure regime where the ion-neutral collision frequency becomes a function of the fluid velocity was described by Godyak [START_REF] Godyak | Soviet Radio Frequency Discharge Research[END_REF] and more recently by Raimbault and Chabert [START_REF] Raimbault | Edge-to-center plasma density ratio in high density plasma sources[END_REF]. Capitalizing on these works, Lee and Lieberman [START_REF] Lee | Global model of Ar, O 2 , Cl 2 , and Ar/O 2 high-density plasma discharges[END_REF] provided a heuristic formula to estimate the one-dimensional h factor in the entire pressure regimes, and this was revisited in Chabert and Braithwaite [START_REF] Chabert | Radio-Frequency Plasmas[END_REF] to provide a simple formula for a 1D Cartesian system of size l: h L = 0.86 3 + 0.5

l λ i + 0.2 T i T e l λ i 2 -0.5 (3)
where λ i is the ion mean free path. This formula has been shown to be in excellent agreement with 1D PIC simulations [START_REF] Lafleur | Edge-to-center density ratios in low-temperature plasmas[END_REF] of inductively coupled plasma (ICP) discharges, for which the electron temperature is fairly independent of space, and will be used as the 1D reference in this paper. It is useful to note that, to our knowledge, there is only one attempt of direct experimental validation of these simplified analytical solutions [START_REF] Kim | Experimental investigation of edge-to-center density ratio in inductively coupled plasmas[END_REF].

When trying to generalize the analysis to two or three dimensions, it may be useful to use a slightly different definition of the h factor, starting from the global particle balance equation:

n 0 n g K iz V = ‹ S Γ i • dS (4)
where n g is the gas density. This equation says that the production of electron-ion pairs by ionization (K iz is the ionization rate) in a plasma volume V has to be balanced by the integral of the ion flux reaching the surface S that encapsulates the plasma. Then, if the integral of the flux is written h n 0 u B S, the h factor becomes:

h = ' S Γ i • dS n 0 u B S (5)
We will use this definition in this paper. We note that in this widely-used definition of the particle balance the following assumptions have been made:

• n 0 in the left-hand side (LHS) of Equation 4is the maximum plasma density in the center of the discharge while it should rather be the average plasma density.

• Electrons and ions are isothermal.

• The dimension of the sheath is small compared to the reactor size.

In this paper we will discuss h factors in 2D ICP discharges with Cartesian coordinates. In section 2, we describe the fluid theory of a non-magnetized plasma discharge for a two-dimensional (2D) Cartesian discharge that will be used as a theoretical support. We have used the work of Kino and Shaw [START_REF] Kino | Two-Dimensional Low-Pressure Discharge Theory[END_REF] as a starting point for our calculation. In section 3, a description of the 2D PIC simulation used to verify the theory is provided. A wide range of pressure is investigated, from 0.3 to 100 mTorr, as well as several aspect ratios. Section 4 compares simulation results with the transport model for electron temperature and h factors. Section 5 explores improvements of the theory and proposes a modified heuristic formula to incorporate 2D effects. Our conclusions are given in the last section.

Theoretical approach of unmagnetized plasma transport in 2D

System definition

The system of interest is a rectangular plasma discharge of dimensions l x and l y , as shown in Figure 1. The plasma is assumed to be sustained by an external power source that does not need to be specified. The system is assumed to be symmetrical with respect to (Ox) and (Oy) axis. In this configuration, from Equation 5, the h factor of the discharge is the following:

h 2D = 2 n 0 u B (l x + l y )    lx/2 0 Γ y (x, l y /2)dx + ly/2 0 Γ x (l x /2, y)dy    (6)
Figure 1: Geometrical configuration of the plasma discharge.

The center of the system of coordinates is taken at the center of the discharge. By symmetry, the fluid velocity is null in the discharge center, and so is the density gradient. The plasma density drops nearly to 0 at the walls, so the simplest boundary condition, known as the Schottky criterion, is the following: n(x, y) = 0 when x = ±l x /2, or y = ±l y /2 [START_REF] Lafleur | Edge-to-center density ratios in low-temperature plasmas[END_REF] This condition is an acceptable approximation to calculate the electron temperature but is not consistent since the flux is finite at the boundary. In addition, it cannot be used directly to evaluate h factors. A better boundary condition is therefore to use the Bohm flux at the boundaries.

General fluid equations

The approach used here aims at generalizing the model investigated in [START_REF] Kino | Two-Dimensional Low-Pressure Discharge Theory[END_REF] to all regimes of pressures, with the assumption that ion-neutral collisions and ionization can be accounted for with constant (i.e. independent of fluid velocity or position) collision frequencies. We will extensively discuss the consequence of this assumption later in the text since it is known to be invalid in the intermediate pressure regime. The momentum conservation equations for ions and electrons are:

m i u∇u = eE -kT i ∇n n -(ν iz + ν i )m i u 0 = -eE -kT e ∇n n (8)
where ion inertia is included and, as in previous models, Boltzmann electrons are used. The continuity equation is:

∇ • (nu) = ν iz n (9)
Expanding these equations on both axis (Ox) and (Oy):

     (u x ∂ x + u y ∂ y )u x = -u 2 B ∂xn n -(ν iz + ν i )u x (u x ∂ x + u y ∂ y )u y = -u 2 B ∂yn n -(ν iz + ν i )u y ∂ x (nu x ) + ∂ y (nu y ) = ν iz n (10)
where it was assumed that T i T e . We look for solutions by separation of variables:

n = n 0 f (x)g(y) ; u = u x (x) u y (y) (11)
where f and g are arbitrary even functions of x and y respectively and n 0 is the plasma density in the center of the discharge. Normalization implies that f (0) = g(0) = 1, and by symmetry, u x (0) = u y (0) = 0. The continuity equation becomes:

u x f f + u x =F (x) + u y g g + u y =G(y) = ν iz (12)
The sum of F (x) and G(y) is constant, which implies that both of these functions are constants which will be called ν x and ν y in the following. Thus, the continuity equation can be separated in two parts:

u x f f + u x = ν x (13)
u y g g + u y = ν y (14)
with

ν x + ν y = ν iz (15)
The momentum conservation equations for each coordinates are:

u x u x = -u 2 B f f -(ν iz + ν i )u x ( 16
)
u y u y = -u 2 B g g -(ν iz + ν i )u y
We can substitute f /f in Equation 16 using f /f = (ν x -u x )/u x , leading to

u x + u 2 B u 2 x (ν x -u x ) + ν iz + ν i = 0 ⇔ u 2 x -u 2 B ν x u 2 B + (ν iz + ν i )u 2 x du x = -dx (17)
such that, after integrating this equation and its equivalent in the y direction, we find:

x = ν iz + ν i + ν x ν 1/2 x (ν i + ν iz ) 3/2 u B arctan ν iz + ν i ν x 1/2 u x u B - u x ν iz + ν i (18)
y = ν iz + ν i + ν y ν 1/2 y (ν i + ν iz ) 3/2 u B arctan ν iz + ν i ν y 1/2 u y u B - u y ν iz + ν i ( 19
)
Since the ions leave the plasma at the Bohm speed and we are in the limit l x , l y s, it follows that:

l x = 2u B ν iz + ν i ν i + ν iz + ν x ν x (ν i + ν iz ) arctan ν iz + ν i ν x 1/2 -1 (20)
l y = 2u B ν iz + ν i ν i + ν iz + ν y ν y (ν i + ν iz ) arctan ν iz + ν i ν y 1/2 -1 (21)
Equations 20 and 21, combined with Equation 15, will set the electron temperature, and the variables ν x and ν y .

As detailed in the Appendix, rearranging the momentum conservation equations and the continuity equations leads after integration to:

f (x) = 1 + u x (x) 2 u 2 B ν iz + ν i ν x -1 2 1+ νx ν iz +ν i (22) g(y) = 1 + u y (y) 2 u 2 B ν iz + ν i ν y -1 2 1+ νy ν iz +ν i (23)
where u x and u y are given by the implicit form of Equations 18 and 19. At the boundaries, x = l x /2 and y = l y /2, the normalized densities in the (Ox) and (Oy) directions -f s = f (l x /2) and g s = g(l y /2) respectively -are given by the Bohm criteria u x (l x /2) = u B and u y (l y /2) = u B :

f s = 1 + ν iz + ν i ν x -1 2 1+ νx ν iz +ν i (24)
g s = 1 + ν iz + ν i ν y -1 2 1+ νy ν iz +ν i (25)
It becomes clear that the edge-to-center plasma density ratio depends on the location where it is evaluated. For example:

h L (y) = n(lx/2,y) n 0 = f s g(y) h L (x) = n(x,ly/2) n 0 = g s f (x) (26)
Using Equation 6and the fact that

Γ y (x, l y /2) = n(x, l y /2)u B Γ x (l x /2, y) = n(l x /2, y)u B ( 27
)
we find:

h 2D = 2 l x + l y    lx/2 0 h L (x)dx + ly/2 0 h L (y)dy    = l x g s β x + l y f s β y l x + l y ( 28
)
where β x and β y are the mean profiles in (Ox) and (Oy) directions respectively. They can be calculated easily by integrating Equations 13 and 14:

β x = 2 l x lx/2 0 f (x)dx = 2u B l x ν x f s ( 29
)
β y = 2 l y ly/2 0 g(y)dy = 2u B l y ν y g s (30)
Finally,

h 2D = β x β y l x l y ν iz 2u B (l x + l y ) (31)
which is consistent with the global particle balance of Equation 4.

Study of the limit cases

2.3.1. High pressure Since in the high pressure case:

ν i ν iz > ν x , ν y (32)
Equation 20 may be simplified since arctan

ν iz +ν i νx 1/2 ≈ π 2 , leading to l x = πu B √ ν x ν i (33)
Using the same approximation along (Oy) gives:

ν x,y = π 2 u 2 B ν i l 2 x,y (34)
such that the condition ν x + ν y = ν iz becomes:

π 2 u 2 B ν i 1 l 2 x + 1 l 2 y = ν iz or π 2 l 2 x + π 2 l 2 y = ν iz D a (35)
Taking this limit in Equations 24 and 25 yields:

f s = ν x ν i = πu B ν i l x g s = ν y ν i = πu B ν i l y (36)
such that the two-dimensional h factor is, from Equation 28:

h 2D = 2u B ν i (l x + l y ) l y l x + l x l y (37)
As shown below, these approximate solutions are consistent with the more classical solutions of the Helmholtz equation that describes the high pressure regime for a 2D Cartesian system [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF]:

n = n 0 cos(k x x) cos(k y y), with k 2 x + k 2 y = ν iz ν i u 2 B ( 38
)
Using Schottky boundary conditions (n = 0 at the plasma edge) and assuming that the sheath thickness s is much smaller than the system size, s l x /2, the following conditions are found:

k x = π l x ; k y = π l y (39)
which yields exactly the same equation as Equation 35for the electron temperature.

Low pressure

At low pressure, the ion-neutral collision frequency is negligible, such that, from Equation 24:

f s = 1 + ν iz ν x -1 2 1+ νx ν iz (40)
For a 1D discharge, l y is infinite, ν y = 0, and ν x = ν iz . Thus, f s = 2 -1 = 0.5, which is the classical result of the low-pressure limit [START_REF] Chabert | Radio-Frequency Plasmas[END_REF], when ionization is included in the momentum balance equation. For a square, ν x = ν y = ν iz /2, and f s = 3 -3/4 ≈ 0.44. This shows that the h factor shows a clear dependance on the aspect ratio of the discharge.

2.3.3.

Large aspect ratios Equation 20 implies that when l x becomes very large ν x becomes very small, such that the following approximation may be used:

l x ν x ≈ πu B ν x ν i + ν iz 1/2 (41)
In addition, taking ν x very small in Equation 24 leads to f s =

νx ν i +ν iz 1/2
. The average of the normalized density profile is then simply β x = 2/π, which is the same as the high pressure case. It is well known that the parameter that sets the shape of the density profile is the product of the pressure by the length. It therefore makes sense that increasing the length in one direction leads to the same effect on β x as increasing the pressure.

General solutions

The analytical formulae derived in section 2.2 can be used to deduce the h factor once the electron temperature has been determined. In this section we explain the procedure used to solve Equations 20, 21, and 15 for the variables T e , ν x and ν y . We introduce the variables θ and φ, defined as

sin 2 φ = ν iz ν iz + ν i (42) tan 2 θ = ν y ν x (43)
so that Equations 20 and 21 take the more compact form:

l x ν i 2u B = cos 2 φG(sin φ cos θ) (44) l y ν i 2u B = cos 2 φG(sin φ sin θ) (45)
with G(x) = (1/x + x) arctan 1 x -1. Arranging those two equations, we obtain:

G(sin θ sin φ) - l y l x G(cos θ sin φ) = 0 (46)
which can be solved numerically to get θ as a function of l y /l x and φ. The relation between the ionization frequency, ν iz = n g K iz , and the electron temperature is gas dependent. In this paper we consider argon, and the ionization rate factor K iz is taken to be [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF]: where T e is expressed in volts. In the first set of results, the ion-neutral collision frequency is taken to be

K iz = 2.
ν i = v i λ i (48)
The gas pressure p, and the gas temperature T g are input parameters of the model. The gas density n g and the ion mean free path λ i are obtained from the following formulae:

p [mTorr] = n g kT g 0.133 = kT g 0.133 σ i λ i ( 49
)
where σ i is the average cross section for ion-neutral collisions, assumed constant and equal to σ i = 1.0 × 10 -18 m 2 , as done in references [START_REF] Lafleur | Edge-to-center density ratios in low-temperature plasmas[END_REF][START_REF] Zhang | Principle of radial transport in low temperature annular plasmas[END_REF]. Note that the spatially constant ion-neutral collision frequency proposed above should be valid at higher pressure (typically p > 100 mTorr for systems of a few centimeters of size), and since ν i becomes negligible at low pressure compared to ν iz , the low pressure regime is also described properly. An improved expression for the ionneutral collision frequency will be proposed in the section 5 in order to better describe the intermediate pressure regime.

At this point, φ is now a given function of the electron temperature, while the pressure and the size of the discharge are parameters. Thus, Equation 44 becomes an equation for the electron temperature which can be solved numerically. This provides solutions for ν x and ν y that can be used in all the analytical expressions of section 2.2. These numerical solutions are compared to PIC simulation results in section 4. Approximate solutions of the model are also possible and detailed in the Appendix, as well as in [START_REF] Lucken | Global models of plasma thrusters: Insights from PIC simulation and fluid theory[END_REF].

2D Particle-In-Cell (PIC) simulations

Description of the LPPic2D code

We have used a 2D parallelized Particle-In-Cell (PIC) code called LPPic2D, developed in-house for the simulation of low temperature plasmas, and described in details in [START_REF] Croes | 2D Particle-In-Cell simulations of the electron drift instability and associated anomalous electron transport in Hall-Effect Thrusters[END_REF]. This code uses a uniform structured mesh, constant time steps, and uses a Monte-Carlo collision (MCC) routine to account for the collisions between the charged species and the neutral gas. Momentum transfer collision cross sections coming from Phelps [START_REF] Phelps | The application of scattering cross sections to ion flux models in discharge sheaths[END_REF] for ion-neutral reactions and Biagi [START_REF] Biagi | Programm magboltz v7.1. cross section compilation. see www.lxcat.net[END_REF] for electron-neutral reactions are implemented. Electron-neutral reactions include elastic collisions, first ionization, and several excitation reactions. Three levels of excitation are resolved for argon. There are two main reactions for ion-neutral collisions which are elastic scattering and backscattering. The code is electrostatic, so only Poisson's equation is solved at each time step, and the magnetic field induced by the motion of the charged particles is neglected. At steady-state, the simulation matches the cell size and time step criteria described in [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF]. Moreover, LPPic2D has been benchmarked using the reference cases published by Turner et al. [START_REF] Turner | Simulation benchmarks for low-pressure plasmas : capacitive discharges[END_REF].

The main assumptions of the model are the following:

• The neutral gas density is assumed to be much larger than those of charged species. Therefore, it is assumed to remain constant and uniform over the simulation runneutral density and temperature are input parameters of the simulation

• Short-range Coulomb collisions between charged species are neglected.

• The dynamics of excited species is not taken into account. Excitation reactions are treated as energy losses.

• Walls are conducting and grounded so they are not charged by impinging particles.

• Finally, secondary electron emission (SEE) is neglected.

A new module was implemented to model an ICP generation. This method was already described by Lafleur and Chabert [START_REF] Lafleur | Edge-to-center density ratios in low-temperature plasmas[END_REF] for 1D simulations of an ICP and extended here in 2D. The heating electric field is uniform over the simulation domain and it is directed along the direction perpendicular to the simulation plane. The absorbed power by the plasma, P abs , is set as an input parameter. We then impose the following electric field,

E = E 0 cos(ωt) e z (50)
and the amplitude E 0 is updated at the end of each radio-frequency (RF) cycle of period T RF using

E 0 = P abs T RF l x l y t+T RF t ˜jz cos(ωt)dxdydt (51)
in order to keep the power absorbed by the plasma constant, where j z is the electron current density along (Oz). The electric field in the simulation plane is derived from the potential, which is the solution of Poisson's equation computed at each time step. The main simulation parameters that were used here are summarized in Table 1.

Figure 2 shows the plasma density obtained in steady-state by the 2D PIC simulation for an argon ICP discharge at 10 mTorr and 19 kW/m 3 . The sheath region, where the electron density drops to nearly zero, while ion density decreases smoothly, is clearly visible in Figure 2(b) and is small compared to the bulk plasma region.

Variable separation

Examining the PIC simulations, it was verified whether or not looking for analytical solutions with separated variables is reasonable. In this prospect, original PIC simulation results for plasma density, n(x, y) were compared with the convoluted density n conv (x, y) = n(x,0)n(0,y) n 0

, where n 0 is the plasma density in the center of the discharge. The model detailed in this paper applies only in the bulk and pre-sheath plasma regions so variable separation in the sheath region is not investigated. The density profiles at the sheath edge -where the discrepancy between the convoluted density and the real one should be the most critical -are plotted in Figure 3 for runs at 1 mTorr, 10 mTorr, and 100 mTorr. The agreement between the real solution and the convoluted one is within the noise (about 5 × 10 14 m -3 ) at 100 mTorr, and it is still very satifsactory at the lower pressures.

Comparison between the fluid theory and the PIC simulations

In this section we compare the electron temperature and the h factor using a constant ion-neutral collision frequency ν i = v i /λ i (Equation 48) in the fluid theory. A large range of pressure is investigated, from the collisionless regime to the higher pressure Schottky regime. We also investigate the effect of the aspect ratio.

Electron temperature

The variation of the electron temperature as a function of l x /λ i is shown in Figure 4. The solid lines are the fluid theory solutions while the dots are the PIC simulation results. Two sets of data are presented. In the first one, l x = l y = 3 cm and the pressure is varied from 0.2 mTorr to 1.0 Torr in order to vary l x /λ i . In the second set of data, the y-direction size is kept constant at l y = 3.0 cm, the pressure is fixed at p = 10 mTorr, Figure 3: Ion density at the sheath edge for various values of pressure. Real profiles are plotted with solid lines and reconstructed profiles with dashed lines. and l x is varied from 0.6 mm to 6 m in order to vary l x /λ i . The curves for the two sets of data cross at l x /λ i = l y /λ i = 10, i.e. for a square discharge. Both curves show a decrease of the electron temperature with l x /λ i , but the decrease is faster when the pressure is varied than when l x is varied.

In both cases a very good agreement is found between the theory and the PIC simulations at low pressure. There is a significant disagreement between the theory and the PIC simulations at the higher pressures, mostly because the electron isothermal assumption breaks down. Indeed, in our PIC simulations, the RF heating electric field is independent of space but since the electron density is not, the power deposition is not uniform. At low pressure, the electron kinetic is largely non-local and consequently the electron temperature remains uniform even if the power absorption is not. At higher pressure, the electron temperature starts to peak in the center, where the absorbed power is the larger, and consequently the assumption of uniform electron temperature is violated.

h factor

The global h factors, labeled h 2D , are computed from PIC simulations and compared to the predictions of the theory described in section 2. In the PIC simulations, the h factors are calculated from the ion current collected at the system boundaries, the central electron density, and the Bohm velocity (computed using the mean electron temperature across the discharge), using the definition provided in Equation 5. Figure 5 shows h 2D as a function of l x /λ i when the pressure is varied, while figure 6 shows h 2D as a function of l x /λ i when l x is varied. The trends are captured by the theory, but the agreement between the theory and the PIC simulation is not very good, in particular when the aspect ratio is varied (Figure 6). The theory curve is also significantly higher than the PIC results in the intermediate pressure regime in Figure 5. As will be discussed in the next section, the theory can be improved by using a more appropriate expression for the ion-neutral collision frequency.

Improved theory and heuristic formula

The intermediate pressure regime

The so-called intermediate pressure regime corresponds to a common situation in industrial plasma applications where the characteristic size of the system is a few times the ion mean free path (1 < l/λ i < 100). Godyak [START_REF] Godyak | Soviet Radio Frequency Discharge Research[END_REF] showed that the ion-neutral collision frequency should, in this case, depend on the local fluid velocity of the ions. This leads to non-linear transport equations, that are difficult to solve analytically, even in 1D. The fluid model used in this paper was already investigated by Sternberg and Godyak [START_REF] Sternberg | Two-dimensional cylindrical plasma for intermediate gas pressures[END_REF] for intermediate pressure discharges in 2D cylindrical coordinates, but no verification of the model was provided, either with experiments or simulations. In the fluid model used so far in this paper, the ion-neutral collision frequency was ν i = v i /λ i (Equation 48) which is a good approximation at high pressure. At very low pressure, the ion-collision frequency goes to zero and therefore the very low-pressure regime (typically l < λ i ) was also well described by the model even when Equation 48 was used. Nevertheless, as seen in Figure 5, the model has to be improved by a better choice of the ion-neutral collision frequency for low and intermediate pressures. We propose to use the form

ν i = v av /λ i (52)
where v av is a typical average velocity of the ions, independent of the position in the discharge, but not independent of the pressure. For all pressure regimes, the ions start at zero fluid velocity in the center and exit the plasma at the Bohm speed. At low pressure, the weakly collisional ions are quickly accelerated and the fluid velocity becomes rapidly higher than the thermal velocity. In this situation, the average ion velocity v av should be close to half the Bohm velocity. At high pressure, the ions drift slowly and the fluid velocity is below the thermal velocity in most of the discharge.

Based on these observations, the following formula was proposed for pressure-dependent average velocity of the ions:

v av = v i + b(α) u B 2 (53)
(α) = tanh 0.2 (α) (54)
This function yielded satisfactory results for argon but another form of function might be more suitable for other types of gases. In the next section, it is shown that the 1D heuristic formula can be straightforwardly extended to the 2D case.

A corrected heuristic formula

The 1D heuristic formula proposed by [START_REF] Lee | Global model of Ar, O 2 , Cl 2 , and Ar/O 2 high-density plasma discharges[END_REF] may be easily corrected to take into account first order 2D effects due to the density profile in the direction perpendicular to the ion velocity vector. As shown in subsection 2.3, this correction factor is 2/π at high pressure. Since the mean normalized density (β x,y ) depends weakly on the pressure and that it has no simple expression, we may keep the value 2/π for all pressure regimes. The corrected heuristic formula then becomes h L,heur,x = 0.55 3 + 0.5

l y λ i + 0.2 T i T e l y λ i 2 -1/2 (55)
h L,heur,y = 0.55 3 + 0.5

l x λ i + 0.2 T i T e l x λ i 2 -1/2 (56)
where 0.55 ≈ (2/π) × 0.86. It should be noted that h L,heur,x -which is the h factor relative to the flux of ions collected on y = ±l y /2 -depends on l y , and not l x . The global heuristic h factor is then h 2D,heur = l x h L,heur,x + l y h L,heur,y l x + l y (57)

Figures 7 and8 show comparisons between the PIC simulations, the two versions of the fluid theory (with different expressions of the ion-neutral collision frequency) and the corrected heuristic formula. The solid dots of Figure 7 are PIC values when varying the pressure between 0.3 mTorr and 100 mTorr. Figure 8 corresponds to several values of l x , from 7.5 mm to 75 mm implemented in the PIC simulation, while l y was always kept to 3 cm. The corrected heuristic formula shows a very good agreement with PIC simulations for both parametric studies because it accounts well for the intermediate pressure regime. The prediction of the original 2D fluid theory described in section 2 agrees reasonably well with the PIC simulation results at low pressure but the discrepancy becomes significant in the intermediate pressure regime. The prediction of the variation of h factor with the aspect ratio is also quite poor with the original theory. However, the agreement is greatly improved when the modified ion-neutral collision frequency (Equations 53 and 54) is introduced in the fluid theory.

Before we conclude this comparison section, we note that at the lowest pressures, it was found in the PIC simulation that the h factor tends to slightly decrease. For these very low pressures, the electron distribution function predicted by PIC simulation becomes strongly anisotropic. This may well be the case in a real system but it could also be an artefact of the heating method implemented in the PIC simulation.

Conclusion and further discussions

The so-called h factors were investigated in 2D Cartesian discharges using a fluid theory and the verification by PIC simulations. The separated roles of pressure and aspect ratio were carefully investigated. The theory predicts the electron temperature with a high degree of accuracy. When the assumption of constant ion-neutral collision frequencies is made, the theory fails in the intermediate pressure regime and the prediction of h factors lacks of accuracy for pressures higher than a few mTorr and below hundreds of mTorr. The theory is greatly improved when a pressure-dependent velocity is used in the ionneutral collision frequency. A corrected heuristic formula, based on the widely used 1D formula, was proposed and it was shown that it matches very well the PIC simulation results in a wide range of pressure and for several aspect ratios of the system. This formula (Equation 55, together with Equations 56, and 57) should be used for global models of a 2D Cartesian discharge. The reasoning developed in this paper applies to cylindrical geometries (see the work of Zhang [START_REF] Zhang | Principle of radial transport in low temperature annular plasmas[END_REF]) and to 3D discharges as well, but with different numerical coefficients.

The PIC simulations have also shown interesting properties on non-ambipolarity. Figure 9 shows electron and ion currents collected on one wall of a square box. The average electron and ion currents are equal, which means that the global flux balance is satisfied, as it must. However, while the ion flux profile follows that of the plasma density, the electron flux features a minimum at the center of the wall and maxima near the sheath edges. PIC simulations carried out here are consistent with a discharge chamber with conductive walls and shows that large currents should be flowing through the walls. Moreover, when the dimensions of the discharge chamber in (Ox) and (Oy) directions are not equal, there are more ions collected on the longer edge and more electrons on the shorter edge, while global current balance is respected. These results are confirming the results found by T. Lafleur [START_REF] Lafleur | Particle-in-cell simulations of ambipolar and nonambipolar diffusion in magnetized plasmas[END_REF]. It highlights the fact that only global h factors should be used in global models, especially for electron power and particle balance. It would be of high interest to conduct these parametric studies again with dielectric walls to see how the profiles of Figure 9 are modified. The heating mechanism could also be improved, for example by solving Maxwell's equations [START_REF] Verboncoeur | An object-oriented electromagnetic PIC code[END_REF], in order to better describe very low pressure regimes with the PIC simulations.

Finally, direct measurements of the h factors remain a challenge but would be of fundamental importance to validate the theory presented in this paper.

was partially funded by CHEOPS project that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 730135. ). The slope of 1.8 is almost not affected when the characteristics of the gas such as the mean cross section for ion-neutral collisions and the ionization reaction rate are changed.

Predicted h factors for one side of the discharge chamber

The advantage of focusing on only one wall of the discharge chamber is to compare the relative role of the pressure and the size of the system in the direction perpendicular to the investigated wall. According to Equation 3, the pre-sheath drop for one side should depend only on the ratio l/λ i . Figure 11, where both parametric studies conducted in this study were plotted together, illustrates that the model developed here (with pressure dependant ion reference velocity, cf Equations 53 and 54) accounts for these distinct trends reasonably well. However, the corrected heuristic formula (Equations 55 and 56) still shows a good agreement with both parametric studies.

Figure 2 :

 2 Figure 2: Electron and ion densities in steady-state for an argon ICP discharge at 10 mTorr and 19 kW/m 3 . (a) Electron density in color plot. (b) Densities of both species represented for x = l x /2.

Figure 4 :

 4 Figure 4: Electron temperature for various pressures and aspect ratios, with l y = 3 cm. The solid lines are the fluid theory solutions while the dots are the PIC simulation results.

Figure 5 :

 5 Figure 5: h 2D factors from the fluid theory (solid line) compared with PIC simulation results (dots) for various values of pressure. l x = l y = 3 cm, p from 0.2 mTorr to 1.0 Torr.

Figure 6 :

 6 Figure 6: h 2D factors from the fluid theory (solid line) compared with PIC simulation results (dots) for various dimensions in (Ox) direction. l y = 3.0 cm, p = 10 mTorr, l x from 0.6 mm to 6 m.

2 i

 2 where b is a function such that b(0) = 0 and lim α→∞ b(α) = 1 which branches heuristically both pressure regimes and α = 4λ lxly is inversely proportional to the pressure squared. As shown in Figures 7, 8 and Figure 11 in the Appendix, good agreements with PIC simulation results were found using a quite steep function: b

Figure 7 :

 7 Figure 7: Predicted h 2D factors compared with PIC simulation results for various values of pressure. l x = l y = 3 cm, p from 0.2 mTorr to 1.0 Torr. The two versions of the fluid model using different assumptions for the ion-neutral collision frequency are plotted together with the corrected heuristic formula (Equations 55 -57).

Figure 8 :

 8 Figure 8: Predicted h 2D factors compared with PIC simulation results for various dimensions in (Ox) direction. l y = 3.0 cm, p = 10 mTorr, l x from 0.6 mm to 6 m. The two versions of the fluid model using different assumptions for the ion-neutral collision frequency are plotted together with the corrected heuristic formula (Equations 55 -57).

Figure 9 :

 9 Figure 9: Electron and ion current densities in steady-state at 10 mTorr and 19 kW/m 3 .

Figure 10 :

 10 Figure 10: Numerical resolution of the model equations with ν i = v i /λ i , for pressures between 0.2 mTorr and 160 mTorr, l x between 0.15 cm and 3 cm, and l y = 3 cm.

Figure 11 :

 11 Figure 11: Predicted h factors for ion collection on the side wall x = l x /2 using the model solution with pressure-dependant reference velocity for the ions (Equations 53 and 54), plotted together with the corrected heuristic formula, and PIC simulation results.

  34 × 10 -14 T 0.59

	e	e -17.44/Te	(47)

Table 1 :

 1 Typical simulation parameters.

	Discharge size in (Ox) direction	l x	0.75 -7.5	cm
	Discharge size in (Oy) direction	l y	3	cm
	Neutral gas pressure	p n	0.3 -100	mTorr
	Neutral gas density	n g	9.6 × 10 18 -3.2 × 10 21	m -3
	Ion mean free path (Ar)	λ i	0.031 -10	cm
	Neutral gas temperature	T g	300	K
	Power absorbed by the plasma	P abs	19.1	kW/m 3
	Excitation frequency	f RF	13.56	MHz
	Number of cells	N cell	40 000 -400 000
	Time step	dt	1 × 10 -11 -3 × 10 -11	s
	Number of CPUs	n proc	200 -400

Acknowledgments

The authors are grateful to A. Tavant and R. Martorelli for the numerous helpful discussions, as well as to A. Lichtenberg and M. Lieberman for their reviews of the manuscript.

This work was granted access to the HPC resources of [TGCC/CINES/IDRIS] under the allocation 2017-A0020510092 made by GENCI, and it

Appendix

Derivation of the model equations

Momentum conservation equation in the (Ox) direction yields:

Introducing the variable s = ux u B √ β , with β = νx ν iz +ν i , and taking into account the fact that s = 0 when x = 0, this leads after integration to:

Replacing with u x , the velocity field is implicitly given by Equation 18.

Since

Hence,

Which leads to the density profile of Equation 22.

Solutions of the model ν x and ν y coefficients represent the importance of ionization in (Ox) and (Oy) directions. We attempted to figure out how the ratio ν x /ν y varies with l x /l y . Numerical solutions of the model were computed for many values of pressure and aspect ratios and plotted in Figure 10. An approximate equation was found under the form:

α is a power coefficient that varies between 1.2 and 2, that depends weakly on the pressure and on the dimensions of the system. In the high pressure limit, Equation 34