Thermal effects in field electron emission from idealized arrangements of independent and interacting micro-protrusions

Modelling studies of thermo-field electron emission (TFE) from protrusions at a cathode surface usually use simulations in 2D axial symmetry. Indeed, time-dependent simulations in 3D are very demanding in computation time. Often, 3D simulations have been restricted to stationary pure field electron emission to account for the drastic current decrease caused by electric field screening when the emitters are close. Little interest has therefore been granted to the heat exchanges occurring between nearby emitters. Although the temperature is a second-order parameter in TFE compared to the electric field, thermal effects become non-negligible in high current density regimes, where self-heating is well established. The present study focuses on the thermal effects occurring during the TFE from micro-protrusions. Our model considers a DC voltage but solves in time the temperature evolution coupling the heat equation and the current continuity equation. The protrusions are modelled as hemiellipsoids with 2D axial symmetry. Emission enhancement due to the increase of the temperature in the thermo-field regime compared to the pure field regime is detailed as a test case for isolated protrusions. Then, full 3D simulations are used to investigate the thermal coupling between multiple neighbouring protrusions via their outwards heat fluxes inside the cathode. The results show a higher current increase due to thermal coupling for dome-like protrusions with a low field enhancement factor. The current increases up to 13% of the total current for aspect ratios of 1, but this value is reached for an extreme applied electric field, hardly reachable in experiments. For sharper protrusions with higher field enhancement, the interaction range through the cathode being shorter, the thermal coupling is suppressed by electrostatic screening. Nevertheless, in arrangements of densely distributed field emitter, when the screening is compensated by a higher voltage, our model predicts the possibility of a moderate but noticeable thermal coupling even for sharp protrusions: a parametric study indicates up to 14.5% of the emitted current being caused by a thermal coupling through the cathode bulk, for protrusions with an aspect ratio of 10 under a fixed applied electric field of 0.4 GV m−1 in DC mode.


Introduction
Field electron emission through quantum tunnelling from a conducting surface, under a strong electrostatic field, is a phenomenon that has led to significant technological developments, enabling the design of a multitude of new electron sources. So-called field electron sources were extensively used in electron microscopy (Swanson and Schwind 2009) and other similar electron beam instruments such as microwave amplifiers (Milne et al 2006) and x-ray generation devices (Sugimoto et al 2010) but also flat panel displays (Spindt et al 1988). Besides, field electron emission is also assumed to be part of the physical process causing vacuum arc ignitions in high-voltage devices that limit, for instance, the performances of particle accelerators (Timko et al 2015). At ambient pressure, the same phenomenon can play a central role in the electrical breakdown of gases at the microscale (Venkattraman 2014, Haase andGo 2016).
Nevertheless, the considerable magnitude required for the electric field to allow field emission, several gigavolts per meter, makes it essentially observable at the apex of sharp conducting protrusions, where the compression of the isopotentials highlights a significantly enhanced electric field (Utsumi 1991, Sarkar andBiswas 2019). Emission from the surface around the apex generates a current inside the protrusion. Above a specific current density, both the Joule heating and the Nottingham effect (Nottingham 1941) can cause a temperature increase, respectively inside the protrusion and at its surface. This self-heating mechanism causes a transition from pure field to thermo-field electron emission (TFE). Higher temperatures then facilitate the electron emission, causing even more ohmic heating. Under a sufficiently high applied electric field, the self-heating mechanism can degenerate into a thermoemissive instability causing a voltage breakdown (Vibrans 1964, Bocharov andEletskii 2007). Otherwise, a steady-state may be reached through an efficient heat sinking towards the cathode bulk, and a cooling Nottingham effect once the Nottingham temperature has been reached (Charbonnier et al 1964). This thermal equilibrium corresponds to steady distributions of current density and temperature.
The phenomenon of electron emission from a fieldenhancing protrusion was already studied in the early 50's (Dyke and Dolan 1956), with building equations initially formulated some decades earlier (Fowler and Nordheim 1928), yet including mathematical errors corrected later (Burgess et al 1953). Due to its complex nature the physics involved could not, however, be fully grasped by analytical calculations and has therefore remained challenging. Two decades ago, numerical simulations have extended the prediction capabilities of the theoretical models established in the twentieth century. Solving systems of coupled differential equations in 1D or 2D with axial symmetry, they have, among others, shed light on the self-heating process of single field emitters (Su et al 1993, Ancona 1995, Rossetti et al 2002. The more recent growth in computing power makes it possible to go further: 3D finite-element solvers allow to track the heating process of multiple emitters and their interactions, up to thermal instability or equilibrium, via time-dependent simulations. For a single emitter, it is clear that the evolution through selfheating at high enough current density from initially pure field emission to thermo-field emission would lead to an increase of emitted current. In the case of multiple emitters in close proximity, the individual emitted current would be damped by electrostatic screening. Under an appropriate electric field, however, high density of emitters could lead in some area to localized heat overflows, which would then enhance the emission via an increase in temperature. This possible coupling of heat fluxes between close protrusions can boost the electron emission and will be here referred to as thermal coupling. The aim of this paper is first to better characterize the transition towards the thermo-field regime over the applied electric field range and then to explore the possibility of thermal coupling. To this end, our simulation routine uses our own TFE model. As in the theory of Murphy and Good (1956), it is assumed that transmission takes place across a planar-imagerounded barrier, often now called a Schottky-Nordheim barrier. The transmission probability through the barrier is evaluated as a function of normal-energy ε n , in a manner generally similar to that used by Murphy and Good. However, instead of their analytical approach towards integration over ε n , the integration is here performed numerically. Doing this enables us to have a single procedure that computes the emission current density at all practical fields and temperatures. It is therefore valid to describe the transition from pure field to thermofield emission regime, while the Murphy and Good system consists of two different formulae for each regime with a gap between them where neither applies. To compute the electric field, our routine uses the Laplace equation. It then solves the time-dependent heat and current equations inside the protrusion. These equations are solved in 3D via a finite-element approach using the COMSOL Multiphysics ® software (5.4) (COMSOL 2020). This way, the routine not only takes into account the electrostatic screening in the inter-electrode space, but also the possible thermal coupling occurring in the cathode volume. Last, the protrusions are modelled by axisymmetric hemiellipsoids. An example of a simulation result for two hemisphere protrusions is shown on figure 1. More details are given in section 2.
The study is divided into two parts, respectively in sections 3 and 4. The first part focuses on single protrusions modelled in 2D: the thermal contributions to the emission is quantified over a broad range of applied electric field in DC, up to the pre-breakdown voltage. This enables to determine the relevant range of parameters where thermal coupling can occur. In the second part, 3D simulations are used to characterize thermal coupling between multiple protrusions, searching for the upper limits to its contribution to the emitted current.

Model
This work takes advantage of the high performance of nowadays numerical solvers of differential equations using a 3D finite-element approach. It is here used to simulate the TFE of well-defined protrusions, distributed at the surface of a cathode and immersed in a homogeneous DC electrostatic field. Example of a 3D simulation performed with COMSOL Multiphysics ® using our thermo-field emission model for two identical hemispheres under an applied electric field of 1.9 GV m −1 . The colormap on the upper surface of the cathode displays the current density variations (colourbar on the left), while the colormap along the vertical cut shows the temperature levels (colorbar on the right).
The corresponding system of equations couples the electron emission with the current continuity equation and the heat balance. It is time-dependent and highly non-linear and therefore needs to be solved numerically in 2D and 3D.

Geometry
The protrusions are here modelled as axisymmetric hemiellipsoids. The main parameters relative to the geometric configuration are thus the height H i and the base radius R i of each ellipsoidal emitter, the corresponding aspect ratio f i = H i /R i , and the positions of their respective symmetry axis in the upper plane of the cathode. The cathode surface, including the surface of the protrusions, is called Σ.
When considering only two emitters, the space of parameters can be reduced to {H 1 , f 1 , H 2 , f 2 , d}, where d is the distance between the symmetry axis of the two emitters (figure 2). Note that R 1 and R 2 are set via f 1 and f 2 respectively.
Protrusions are placed in a volume domain bounded by ±L x in the x-direction, ±L y in the y-direction, and {−L cat , L z } in the z-direction. These bounds are chosen so that they do not influence the physics, but their impact on computation time is limited. A good compromise is obtained with L x , L y , L z and L cat of the order of a few dozen times the height of the highest protrusion.

Details of the electron emission model
• Following Sommerfeld theory of metal, the electrons inside the protrusion are described by a Fermi-Dirac statistic.
Their energy distribution depends on the protrusion equilibrium temperature T and the material work function φ. • In the inter-electrode space, the solving of the electrostatics gives the electric field distribution over the protrusion surface Σ. • The corresponding electric field magnitude F is taken into account together with the image charge correction to compute the Schottky-Nordheim potential barrier. • The transmission probability D(F, ε n ) of an electron of normal-energy ε n through the barrier is computed by solving the Schrödinger equation in 1D using the Kemble formalism (Kemble 1935). Curvature effects are therefore not taken into account, as they have been shown to be negligible for apex radius larger than a few nanometers (Kyritsakis and Xanthakis 2016). • The mathematical formulae obtained are then computed numerically to get the emitted current density J e (F, T, φ) from D(F, ε n ). More details are given in appendix A. If the electric field is high enough over the protrusion apexes, electrons around the Fermi level are able to escape the metal through quantum tunnelling, and a current is extracted. • The emitted current causes a Nottingham effect at the protrusion surface. The corresponding Nottingham heat flux Φ N is computed by accessing the mean energy difference W N between the emitted electrons in vacuum and the replacement electrons from the bulk 3 : • At the same time, inside the protrusions, the establishment of a current leads to Joule heating. As a consequence, the temperature increases along the protrusion axes. If the initial current density is high enough, the emitter may then enter the thermo-field emission regime where the increase in temperature brings many electrons towards higher normalenergies, facilitating the emission and yielding even more extracted current. • This positive feedback loop between Joule heating and current is then counterbalanced over time by a heat sinking towards the cathode bulk-assumed here acting as a thermostat kept at 300 K-and in addition a cooling Nottingham effect if the Nottingham temperature is exceeded. • When a steady state is reached, the equilibrium distributions of the local variables, such as the current density and the temperature, are accessible as well as the global variables like the intensity of the total emitted current.

Electrostatics model
The electrostatic field at the emitters surface is computed solving the Laplace equation for the potential V out in the interelectrode gap, assumed to be perfect vacuum: where V out and ⃗ F respectively denote the electric potential and field. The corresponding boundary conditions (BCs) are the following. At the cathode surface Σ (see the green dashed line in figure 3), a Dirichlet BC sets the V out to the macroscopic applied potential V 0 . At the upper limit of the simulation domain for z = L z a von Neumann BC sets the field norm F to the applied electric field E 0 = V 0 /D gap , where D gap is the distance separating the anode from the cathode. On the far left and right sides, where x = ± L x or y = ± L y , symmetry von Neumann BCs are enforced on the electric field.
From the local electric field norm F at position ⃗ r on the protrusion surface, the field enhancement factor (FEF) is defined by: while its specific value at the apex of a protrusion (apex FEF) writes: with F a the field value at the apex. The FEF is defined for an isolated protrusion, i.e. not interacting with any other protrusions. It is defined here as an intrinsic property of the protrusion, directly related to its geometry.
Note that field modification due to space charge during the emission is neglected at present. The interested reader can consult (Dyke et al 1953) (experimental) or (Zhu and Ang 2015, Seznec et al 2020) (modelling).

Electric model
Inside the cathode, charge conservation must be satisfied: with the current ⃗ j expressed in the form of the gradient of V in : where V in denotes the local potential inside the cathode and its gradient is obviously related to the electron emission. The equation resulting from (6) and (7) is solved with the following BCs: Note that according to (7), the electron emission J e at Σ implies a small voltage loss along the axis of the protrusions so that V in | apex ⩽ V in | Lcat . This loss is amplified by higher current density and electric resistivity as often observed for carbon emitters (Minoux et al 2005). In our case, however, simulations use typical metal conductivities of the order of 10 6 -10 7 S m −1 . Moreover, as a distance of D gap = 200 µm is chosen between the two electrodes, the voltage applied in our simulations is around 10-100 kV. The relative loss of voltage is therefore completely negligible, of the order of 10 −6 , leading to V in | apex ≲ V 0 . This is why the electrostatics inside the enclosure can be solved only once assuming V out | Σ = V 0 all over the cathode surface, as presented in section 2.3.

Heat transfer model
The electron emission causes Joule heating inside the protrusion and Nottingham heat fluxes at their surface. The thermalization of the cathode is then described by the following heat equation: with T the temperature, ⃗ ϕ the corresponding heat flux, µ, c, κ and σ respectively the volume mass density, the specific heat capacity, the thermal conductivity and the electrical conductivity of the cathode material, all being temperature-dependent.
In all the results of the present paper, the cathode bulk is assumed big enough to act as a thermostat beyond L cat . At the interface Σ, cooling or heating flux condition is fixed by the Nottingham effect Φ N . On the far left and right sides of the simulation domain, symmetry von Neumann BCs are set. These BCs are summarized by the following equations: Note that at the emitter surface, no radiation losses are taken into account as they were found to be negligible compared to the conductive heat evacuation at the protrusions bases. Over the parameters range of interest, an upper estimation can be obtained assuming black body radiation: where σ ⋆ = 2π 5 k 4 15h 3 c 2 is the Stefan-Boltzman constant.

Simulation algorithm
A complete simulation corresponds to the chart-flow depicted in figure 4. To get the electric field at the protrusions surface, the electrostatics are only solved once using a static finiteelement method (FEM) included in COMSOL Multiphysics ® (5.4) (COMSOL 2020). With this software, the coupled electric and heat equations are then solved together using a timedependent FEM and a non-linear solver (Newton-Raphson) to get the current density and temperature distributions. These distributions are updated at each time step by computing the value of the emitted current density J e (F, T, φ) and the Nottingham heat flux Φ N (F, T, φ) for the BCs, via our electron emission model written in Fortran. The simulation stops when a steady-state is reached. The model is versatile in the sense that it is able to follow the electron emission of any 3D configuration of any number N of protrusions of any shape, of course at the cost of a computation time increasing with both N and the shape complexity.

Thermal effects in the field emission of a single protrusion
The results presented here focus on the influence of thermal effects in field electron emission from idealized protrusion with hemiellipsoidal shape. Comprehensive insights on the heating process of a single protrusion are given, in order to characterize the transition from the pure field emission regime to the thermo-field emission (TFE) regime. Because only one emitter is here considered, the x and y coordinates described in the model section are merged in a radial coordinate r and the simulation domain is reduced to a 2D box with axial symmetry (see the model described in Seznec et al 2016, section 2 therein). This dimension reduction eases the exploration of thermal effects over an extensive range of applied electric field. Besides, it offers the possibility to extract the more relevant parameters to study then thermal coupling (section 4), which on the contrary do require a 3D modelling.
In the present work, the space of parameters is reduced to the aspect ratio f by setting the protrusion height to H = 10 µm. The results focus on aspect ratio from 1 to 10, especially f = 1, f = 2, f = 5, and f = 10 in this section, corresponding respectively to an apex FEF of β a = 3, β a = 5.76, β a = 17.9, and β a = 49.3. Besides, our simulations are usually performed with three metals, namely tungsten, molybdenum, or titanium. Hereafter in this document, all the results presented are obtained with titanium: average work function φ = 4.3 eV, melting temperature T M = 1941 K, and conductivities σ(T) and κ(T) set according to the fitting polynomials proposed for titanium in Milosević and Maglić (2008) and Ho et al (1972) respectively.

Transition to the thermo-field emission regime
The details of the variation of the current with the voltage is difficult to predict and a complete simulation is therefore usually required for each different arrangement. The causal link between the emission and the voltage, however, is simple: the higher the applied voltage V 0 is, the more current is extracted. In particular, one can reasonably predict that at some applied voltage, the current density inside the protrusion should be high enough to induce a noticeable increase of temperature via the Joule heating and/or a heating Nottingham effect. There would then be a transition over the applied voltage from pure field electron emission (low current density, constant temperature) to thermo-field electron emission (high current density, self-heating process). If the voltage-and so the related current density-is high enough, various thermal instabilities may develop as a result of the self-heating process that may degenerate into vacuum breakdown. For instance, as the temperature at the protrusion surface rises, field-enhanced surfaceatom diffusion may become significant and cause the growth of additional nano-tips at the protrusion surface. According to Kyritsakis et al (2018), such a phenomenon could initiate a thermal runaway leading to partial evaporation of the protrusion. However, except for refractory metal (Jansson et al 2020), the temperature required to cause a significant surfaceatom diffusion generally coincides with the melting temperature. As our model does not include the physics that would describe such evolutions towards a vacuum breakdown, it is only valid up until the pre-breakdown condition is met. Therefore, hereafter in this work, the pre-breakdown condition is assimilated to the maximum temperature reaching the melting temperature, anywhere in the protrusion. This occurs for a threshold electric field denoted E pb 0 that is produced by an applied voltage considered hereafter-for simplification-as the pre-breakdown voltage: Thus, E pb 0 is the minimum applied field for which the melting temperature is reached somewhere in the cathode (T M = 1941 K for Titanium). In order to guarantee the integrity of the protrusions, the field values in the simulations are restricted to E 0 ⩽ E pb 0 : the applied voltage never goes beyond the pre-breakdown voltage.
The emitted current I pb when E pb 0 is applied is therefore the maximum current one can extract from a single protrusion before it begins to melt. This field value is therefore used to define the emission range of a single protrusion is defined as an arbitrary lower limit for which the protrusion emits one thousandth of I pb .
Over this emission range, the transition from the cold to the thermo-field regime is explored for four typical protrusions, f = 1, 2, 5 and 10, via the following global quantities: • The self-heating of the protrusion in terms of the mean value of the Joule and Nottingham heating in relation to the protrusion volume V: • The maximum temperature at equilibrium in the protrusion volume: • The total emitted current: • And the relative gain in current due to the increase in temperature: All these global quantities introduced by equations (12)-(15) are depicted versus the applied electric field in figure 5. Taking advantage of the different aspect ratios considered, the results for each protrusion are displayed on the same x-axis. Each data point at a given E 0 corresponds to a 2D axisymmetric transient simulation whose quantities are evaluated once a steady-state has been reached. Note that for the parameter range presented in this work, the convergence time towards a steady state is between a few microseconds to several dozens of microseconds, which is in accordance with the rough estimate τ ∼ H 2 /α that one get from the material thermal diffusivity α = κ/µc.
The first panel of figure 5 shows the respective contribution of the Nottingham effect and the Joule heating in the thermal equilibrium at different applied electric fields. It is interesting to notice that for low and medium values of applied field, the Nottingham effect contributes alone to the heating. This is consistent with the results of Su et al (1993) 4 and Ancona (1995) 5 . It is only at high electric fields (i.e. high enough emitted currents) that the averaged Joule heating catches up, where the averaged Nottingham heating reaches a maximum. For even higher field values, the averaged Nottingham effect decreases (and can then be cooling) while the ohmic heating takes over and bring the tip closer to breakdown. Note however that the overall heat surplus P J + P N is always compensated by the for aspect ratios f = 1, 2, 5 and 10.
• The first panel gives the variation of the averaged Joule heat P J and Nottingham heat P N . The total heat P J + P N is explicitly displayed at g I = 1% (see bottom panel, hatch separation) and at E pb 0 . A vertical dashed line highlights the Nottingham inversion. Note that the y-axis is in symlog scale: linear between −0.01 and 0.01 and logarithmic beyond. The insets show the shape of each considered protrusion with its aspect ratio, apex FEF and volume.
• The second panel shows the increase of the maximum temperature Tmax inside the protrusion. The temperatures at g I = 1% and at the Nottingham inversion are explicitly displayed. • The third panel presents the current emitted from the protrusion assuming cold emission only: I(F, 300 K) and the total current when the thermo-field emission supplement the current: I (F, T). The values at E min 0 and E pb 0 , as well as the value at g I = 1% are explicitly displayed. • The last panel gives the relative gain in current due to the temperature increase, i.e. the relative discrepancy between the above curves for I(F, 300 K) and I(F, T). A separation at g I = 1% is highlighted by blue and red hatches to arbitrary delimit the pure field regime from the thermo-field regime. The corresponding field value is given in line with the pre-breakdown field E br 0 for each aspect ratio f.
heat sinking, as all of the presented results come from simulations that have reached a steady-state. The second panel shows the continuous increase of the protrusion maximum temperature towards the melting temperature T M = 1941 K, as a result of the self-heating process.
The field values E Nh 0 for which P N = 0 are highlighted by extra x ticks on the first panel and the corresponding maximum temperatures T Nh are given on the second panel (see the dashed gray lines). The values are also summarized in table 1. Remarkably, the averaged Nottingham effect reverses at smaller temperature for smaller aspect ratio f. This is the consequence of a smaller local electric field at the emission surface as will be shown in section 3.2.
To assess the contribution to the current of the increase in temperature, two sets of results are compared in the third panel. The first set corresponds to electrostatics-only simulations, i.e. cold emission only. In this case the current I(F, 300 K) is obtained by setting the temperature at 300 K all over the emitting surface (blue crosses-figure 5-3rd panel). On the contrary, the second set comes from a complete solving of the model, as described in figure 4, and the computation of the current I(F, T) uses the whole equilibrium temperature distribution at the emitting surface (red circles-figure 5-3rd panel). When the applied electric field is much lower than the pre-breakdown field E pb 0 no distinction can be made between the two curves. However, closer to E pb 0 a difference appears that witnesses a gain in current g I due to the increasing contribution of 'hot' electrons that cross the barrier at normalenergies significantly above the Fermi level. This corresponds to the thermo-field emission regime.
The last panel displays g I (E 0 ) which is the relative discrepancy between the red curve I(F, T) and the blue one I(F, 300 K) of the third panel. For f = 10 the gain on the current goes up to g I = +61% over the applied electric field range ∆E 0 = [0.095, 0.16] GV m −1 . For f = 1 it reaches g I = +109% over a relatively smaller range, ∆E 0 = [1.3, 1.9] GV m −1 . This tends to show that the rounder the protrusions are, the more thermo-field effects contribute to the emission.
Besides, g I (E 0 ) is the perfect quantity to distinguish the pure field electron emission regime (g I ∼ 0) from the thermofield regime (g I ⩾ 0), as it directly quantifies the contribution of the thermal effects to the current. By setting an arbitrary limit between the two regimes at g I = 1%, a separation appears for each f at applied field between 0.75 and 0.80E pb 0 with a corresponding heating around ∼5 µW µm −3 . The applied field corresponding to this separation is denoted E 1% 0 . Because the self heating of protrusions is more directly related to the current-both for Nottingham and Joule heating-it is interesting to read this limit in terms of current per unit volume. Dividing the current values given on the third panel at the 1% separation by the protrusions volume V shown on the insets of the first panel, one gets a similar value for each protrusion, between 25 and 30 µA µm −3 (thermal effects are still too small to make a clear difference). Doing the same computation at the pre-breakdown voltage, however, leads to a more readable distinction between protrusions, I pb /V increasing from 810 µA µm −3 for f = 10 to 1700 µA µm −3 for f = 1. Values are given in table 1. From this result, it appears that rounded protrusions require more current per unit volume to heat up to the breakdown. This can be mainly related to the fact that rounder protrusions have a better heat sinking at their base, also meaning a temperature increase that radially spreads wider. This assertion is developed in the next subsection.

Details of the emission in the thermo-field regime
The self-heating process has been shown to differ between protrusions depending on their aspect ratio f. To better understand this, it is of interest to look at the distributions of the local variables F, T, and J e (F, T) at E 0 = E pb 0 . Figures 6 and 7 displays these distributions for f = 1, 2 and 5.
The first panel of figure 6 shows the field distribution at prebreakdown voltage, so for an applied field E 0 = E pb 0 . Although E pb 0 is higher for lower aspect ratio, the graphs highlight that the local field value at the apex F a is smaller for lower aspect ratio while the field enhancement spreads over a larger surface. Hence, for smaller aspect ratios, the emission occurs at lower field values but over a larger surface around the apex. This is measurable via the introduction of the radius R 90% defining, for an axisymmetric situation, the protrusion surface around the apex that emits 90% of the total current: being the total extracted current and  2D axisymmetric emission of single protrusions with aspect ratio f = 5, f = 2 and f = 1 when the corresponding pre-breakdown field E pb 0 is applied. The upper schemes recall the protrusion properties. The graphs of the first row compare the radial distribution of the local electric field over the whole protrusion, from r = 0 to r = R with a common y-axis. Fa corresponds to the local electric field at the apex. The second row of graphs presents the distribution of the emitted current density, taking into account the self-heating process (thermo-field regime: Je(F, T)) or not (pure field regime: Je(F, 300 K)). Once again, the y-axis is the same, but for better readability, current densities have been scaled by a factor of 4 and 10 for f = 2 and f = 1 respectively, as indicated in bold in the legend. Note that the distributions are zoomed (dashed grey lines) on the area emitting 90% of the total current. These areas, translated into radial coordinates, are highlighted by red and blue hatches. Isothermal curves for the 2D axi-symmetric emission of single protrusions with aspect ratio f = 5, f = 2 and f = 1 when the corresponding pre-breakdown field E pb 0 is applied. The step between two successive isothermal curves is 100 K. The thick dotted lines above the protrusion apexes highlight the area contributing to 90% of the emitted current. R 90% is the corresponding radius.
The overall lower electric field over the emission surface for smaller aspect ratios explains why rounder protrusions benefit more from thermo-field contribution (higher g I at E pb 0 ). Indeed, the influence of temperature on J e (F, T) predicted by the theory is precisely higher at lower fields.
The larger emission area for smaller aspect ratio explains why their emitted current at pre-breakdown voltage is higher while their current densities are lower: approximately four times lower for f = 2 compared to f = 5, and ten times lower for f = 1, as highlighted in the legends of the current density graphs of figure 6. The emission area being larger for smaller aspect ratios, the Joule heating takes place over a larger volume beneath it. The increase of the volume temperature thus extends wider and can reach the surrounding of the protrusion basis. This is particularly visible for f = 1 and f = 2 on figure 7.
The question then arises as to whether this heat diffusion for low aspect ratio protrusions could significantly impact the emission of other protrusions in the vicinity. Characteristics of this possible thermal coupling are studied in the following section.

Phenomenology of interactions in 3D
So far, simulation results of single protrusions have been analysed. In real situations, as observed in electron sources or with surface roughness, the emitting protrusions are usually grouped. In such cases, even though each emitter still might be assumed to have axial symmetry, the whole situation is not symmetric anymore. 3D modellings of the emission are required to account for the proximity of the protrusions. In the present work, two simplified arrangements are considered: • A two identical protrusions arrangement separated by a distance d, named hereafter the 2IP arrangement ( figure 8(a)). 2IP stands for Two Identical Protrusions. • An infinite 2D lattice of identical protrusions with a spacing distance d, named hereafter the array arrangement ( figure 8(b)). The unit cell of the array is a d × d square.
The whole arrangement is however modelled from a d 2 × d 2 square including a quarter of a protrusion, and using mirror symmetries at the lateral boundaries ( figure 9).
In such configurations, the proximity of the protrusions induces electrostatic screening. This simple, well-documented phenomenon (Dall'Agnol et al 2018)-also named mutual depolarization (Forbes 2016)-reduces overall the field enhancement capability of the protrusions and consequently the local electric field around the apex. This effect is here quantified by the relative loss of the apex electric field δ a : where F a is the actual electric field at the apex of the interacting protrusion (i.e. with screening) and β a E 0 corresponds to the apex electric field for a single isolated protrusion (i.e. without screening). Due to the exponential decrease of the emitted current with smaller electric fields, even a small field loss induces a noticeable current loss. This is particularly troublesome for electron sources, leading to many studies about the spacing optimization in field emitters array (FEA) (Dionne et al 2008, Bieker et al 2019, de Assis et al 2020. However, these studies often do not consider thermal effects in their search for an optimal spacing. In the thermo-field regime, electrostatic screening reduces the emitted current and therefore reduces the self-heating process as well. As an example, figure 10 compares for f = 1 the temperature map of a single protrusion ( figure 10(a)) and this of a 2IP arrangement with d = 3R ( figure 10(b)) under the same applied electric field: E 0 = 1.9 GV m −1 . The decrease of the maximum temperature from 1941 to 1471 K reveals a significant self-heating reduction for a field loss at the apex of δ a = −2.5% only. As a consequence of this electrostatic screening, close protrusions need a higher electric field to reach the thermo-field regime.
It is however remarkable that a second effect, also related to proximity, slightly compensates the temperature decrease caused by screening. This is what is here referred to as thermal coupling and is highlighted on figure 10(b) by the joining of the isothermal curves in the cathode below the protrusions. An indirect way to evaluate its influence on temperatureand so on current-is to systematically compare two sets of independent simulation results. A first set corresponds to a standard simulation in 3D, accounting for both electrostatic screening and thermal coupling. For the second set, the electrostatics are first solved in 3D to account for the screening, but the heating process is then treated independently for each protrusion 6 . Proceeding this way artificially cancels the thermal coupling. The corresponding results are thus labelled thermal coupling off (TC-off) while the results of the first set are labelled thermal coupling on (TC-on). The consequence of cancelling the thermal coupling for the previous example is shown in figure 10(c). It highlights a loss of 35 K in the maximum temperature and about 100 K in the temperature at the protrusion base compared to figure 10(b), highlighted by the position of the isothermal curve at 700 K (displayed in white). Although this is a rather small variation, the thermal coupling is clearly shown in this example.
In fact, the small magnitude of thermal coupling is due to the electrostatic screening. To support this statement, a third simulation is performed. It independently solves the electrostatics for each protrusion alone and then uses the resulting surface electric field to simulate the heating process once including all protrusions. This procedure artificially cancels the electrostatic screening and the corresponding results are labelled electrostatic screening off (SC-off). The consequence is shown in figure 10(d) which is to be compared with the single protrusion case ( figure 10(a)). The effects of thermal  coupling without screening are here roughly twice higher: the maximum temperature increases by 70 K and the increase at the protrusion bases is rather around 200 K.
This example gives insights into the phenomenology of thermal coupling. More details on the parameter dependence of thermal coupling are given in the following.

Parametric analysis of thermal coupling
The goal of this subsection is to explore the parameters that influence thermal coupling and find the range that maximises the effect on current. A priori: rounded protrusions, close to breakdown, densely distributed so that they heat each other but spaced enough so that the screening does not suppress the emission. Because this work aims at finding an upper limit to this phenomenon, analyses are first performed on protrusions with an aspect ratio of 1, as dome-like protrusions were shown in section 3 to be the most favourable case for the heating of the cathode and the contribution of temperature to the current. This is also the reason why the array arrangement was chosen-in addition to the 2IP arrangement which is the easiest arrangement of multiple protrusions. Nevertheless, f = 1 protrusions require a particularly high applied electric field to emit in the thermo-field regime (E 0 ⩾ 1 GV m −1 ). The following analysis should therefore be considered as a theoretical investigation of the thermal coupling phenomenon, in the range of parameters that maximizes it. Section 4.3 will be dedicated to results and discussion in more practical situations.
To get a quantitative evaluation of thermal coupling influence on the emitted current, both TC-on and TC-off sets are compared: the discrepancy between the current per protrusion in both sets, respectively I on and I off , is denoted ε I and writes: Because I on is taken as the reference, ε I is a measure of the contribution of thermal coupling to the current. For instance, ε I = 4% means 4% of the current I on is due to thermal coupling. In order to get an upper limit, the aim is to find for both the 2IP and the array arrangement the distance d that maximizes ε I . The results are compared in figure 11 for an applied electric field E 0 = E pb 0 (f = 1) = 1.9 GV m −1 , i.e. an apex electric field F a = β a E 0 = 5.6 GV m −1 .
Firstly, it is essential to note the significant difference in screening between these two arrangements. The screening effect appears to be much higher for the array than for the 2IP arrangement. This is understandable as an infinite array becomes quite close to a flat surface when all the tips are in contact. For the array arrangement, when contact between the protrusions is reached, on the one hand, the maximum field loss is δ a = −33%. The emitted current consequently drops from 3.55 A to 2 mA only. No temperature increase is produced in this case, and thermal coupling is null. Beyond 16R on the other hand, no more field loss nor thermal coupling is possible and the system behaves like a sum of isolated protrusions. Therefore, the spacing d maximizing the thermal coupling is reached for a compromise in the distance, here d = 5R, still allowing the heat fluxes to couple while preventing a too high screening: δ a = 5.2% allowing a current of 1.3 A, in the thermo-field regime, of which 13% are due to thermal coupling. For the 2IP arrangement, the maximum field loss does not exceed −4.7%. The maximum coupling can thus be reached when the two protrusions are in contact. The emitted current I = 1.2 A is similar to the array, but thermal coupling contributes here for 3.5% only. A higher impact of thermal coupling was found for the array arrangement, as expected: the more crowded the protrusions, the smaller the dissipation volume.
Overall, however, the observed values for ε max I are rather small in both cases since, having been obtained for f = 1, they are assumed to be upper limits. This confirms that when dealing with thermo-field emission, the thermal coupling due to protrusions proximity only has a second-order effect on the emitted electrons. The electrostatic screening plays the major role and preventssignificant thermal coupling. In fact, by artificially cancelling the screening in the simulations, all else being equal (results not shown), the contribution of thermal coupling gets roughly doubled: ε I at 2R = 7.9% for the 2IP arrangement and ε I at 5R = 25% for the array arrangement (not being the maximum anymore).
To go further, one can search for the maximum current discrepancies ε max I over d for higher f and analyse their behaviour. The results are limited to the array arrangement as the thermal coupling for a more sparse arrangement-all else being equal-would be smaller. As for f = 1 (see the graphs of figure 11), pairs of simulations with (TC-on) and without (TC-off) thermal coupling are performed for higher f over a Figure 11. Comparison of the influence of thermal coupling on the emitted current per protrusion between the 2IP and the array arrangement for f = 1. The x-axis corresponds to spacing distances d ranging from 2 to 16 protrusions radii and is the same for all graphs. The first row presents the screening efficiency at the apex δa. The middle row compares the emitted current with and without thermal coupling and the bottom row exhibits the discrepancy between those two cases. This discrepancy is in fact the percentage of emitted current that is due to the thermal coupling between the protrusions (ε I ). The maximum discrepancy over the tested distance d is highlighted by red hatches, and is denoted ε max I .
wide range of array spacing d, with a sampling step equal to the protrusion radius R. The applied electric field corresponds respectively to the pre-breakdown field of each aspect ratio found in section 3 (figure 5). The value for ε max I is obtained from the pair exhibiting the highest current discrepancy, and the corresponding spacing at maximum is denoted d * . Results are reported in figure 12, and in table 2 along with the field values, spacings at maximum, and maximum temperatures. It is observed that the thermal coupling contribution to the current quickly becomes negligible for higher aspect-ratio: from a maximum contribution over d of 13% for f = 1, ε max I drops to 2.2% for f = 4. In fact, as observed in figure 7, sharper protrusions dissipate their heat over a smaller distance at their bases. A thermal coupling can therefore only occur for shorter array spacing, where the higher electrostatic screening reduces the thermo-field emission. This assertion is supported by the spacing at maximum that decreases from 50 µm for f = 1 to 25 µm for f = 4, almost doubling the apex field loss: from −5.23% to −10.6% (table 2).
For f ⩾ 4 (results not shown), the thermal coupling becomes too weak to be observed at large distances (⩾30 µm), and is suppressed by electrostatic screening at small distances (⩽25 µm). In between, several local maxima appear as distinct trade-offs between a higher proximity or a smaller screening. Overall, these maxima correspond to negligible contribution to the current, ϵ I ⩽ 2%, and disappear for f ⩾ 8.
These results highlight the difficulty to observe a thermal coupling between sharp protrusions. However, the applied voltage was here limited to the pre-breakdown voltage of isolated protrusions. With neighbouring protrusions in interaction, because of electrostatic screening, the system in fact requires a higher applied electric field for the protrusions to reach an equivalent thermal situation. The actual maximum temperatures for the data of figure 12 are below the melting  temperature: T on max = 1316 K for f = 1 and T on max = 962 K only for f = 4 (table 2). By ramping up the voltage to its effective pre-breakdown value, the values of ε I for high aspect ratio could then become noticeable. It is therefore hard to set an upper limit on thermal coupling magnitude without setting an upper limit on the applied electric field. This issue is discussed in what follows.

About the magnitude of thermal coupling in experiment-like situations
In modelling, physical processes can be isolated from one another. In addition, there is no intrinsic limit to the magnitude of the parameters. In experiments, however, the magnitude of the applied electric field usually ranges well below the gigavolt per meter. On the one hand, indeed, the maximum voltage delivered by the generators is limited. On the other hand, a wide variety of physical processes can couple and cause breakdown at lower field values than the ones predicted only by the vaporization of an idealized micro-metric protrusion. In fact, in experiments concerning field electron sources, the applied field usually ranges from several volts per micrometers for the modern highly enhancing structures such as carbon nanotubes (E 0 ≳ 10 6 V m −1 ), to a few tenths of gigavolt per meter for less sharp structures (E 0 ≳ 10 8 V m −1 ) (Li et al 2015). In the case of high-gradient accelerating structures, modern linear accelerators such as CLIC-the Compact Linear Collider at CERN-operate with RF fields, around the 100 MV m −1 range (Wang et al 2012). Several experimental studies evaluate the high-voltage resistance of accelerators in DC mode and exhibit breakdown fields from 100 to 800 MV m −1 (Descoeudres et al 2009).
With these values in mind, here is finally investigated the magnitude of the thermal coupling under a purposefully borderline field of 0.4 GV m −1 . Under this field restriction, the upper limit of ε I over our parameter range is estimated. Referring to the third panel of figure 5 in section 3, extracting current via thermo-field emission at 0.4 GV m −1 from protrusions with f ⩽ 5, i.e., β a ⩽ 17.9 is not possible. Hemiellipsoids from which a current can be extracted have an aspect ratio f ⩾ 5. Therefore, a final sequence of simulations is performed for f = 5-10 with E 0 fixed to 0.4 GV m −1 . Note that this value is much higher than the pre-breakdown voltage of isolated protrusions of aspect ratio f ⩾ 5. In an array arrangement, however, because of the electrostatic screening, there exists for each of these aspect ratios a unique spacing distance for which the maximum temperature in the protrusions reaches the melting temperature. In analogy to our definition of the pre-breakdown voltage V pb , this pre-breakdown spacing is denoted d pb . For this specific spacing, two sets of simulations in the array arrangement are compared following the method described previously: the first set takes into account the thermal coupling while this phenomenon is artificially cancelled for the second set. The relevant results for the comparison are summarized in table 3.
Obviously, because sharper protrusions have a higher apex FEF (denoted as β a ), they can emit in the thermo-field regime while undergoing a higher screening than low aspect ratio protrusions. This is why the distance d pb is smaller for higher f. For example, for f = 6 the pre-breakdown is reached at d pb = 16.1 µm with a screening of δ a = −21.5%, while for f = 10, d pb = 5.8 µm is smaller, inducing a roughly three times higher apex field loss, δ a = −60.4%. Despite the smaller impact of thermal effects for sharper protrusions exposed in section 3, the applied electric field of 0.4 GV m −1 allows higher proximity for higher f, enhancing further the influence of thermal coupling for f = 10: ε I = 14.5% than for f = 6: ε I = 6.2%. This remarkable result confirms that the main obstacle to a noticeable contribution of thermal coupling to thermo-field emission is the electrostatic screening. When the latter phenomenon is counterbalanced by a high enough applied field, the contribution of thermal coupling can be noticeable. In the present case, ε I ≳ 10% with d ≲ H for sharp protrusions with apex FEF corresponding to regular field emitters characterized by 20 ⩽ β a ⩽ 50. Such configurations can be observed in FEAs experiments where the emitters are grown densely packed. Table 3. Comparison at the pre-breakdown distance d pb of the maximum temperature and emitted current between a set that takes into account the thermal coupling (T on max , Ion) and a set that artificially cancels thermal coupling (T off max , I off ). The corresponding contribution of thermal coupling to the current (ε I ) is given. The apex field Fa and apex field loss δa are displayed as well for more clarity. The comparison is limited to protrusions with aspect ratio f from 10 to 5, as no thermo-field emission is possible for rounder protrusions (f ⩽ 5) under an applied electric field E 0 = 0.4 GV m −1 . The protrusions height is 10 µm. For example, in Spindt (1992) the authors report two different FEA made of molybdenum and silicon dioxide with emitter height H ≳ 1 µm and H ≳ 10 µm, respectively separated by a distance d ∼ 2 µm and d ≳ 10 µm. It is possible to roughly evaluate the effective FEF to 10 1 ⩽ β ⩽ 10 2 from the experimental I-V characteristics. Besides, the order of magnitude given on the emitter curvature radius makes it possible to approximate the range of emission current density per emitter: 10 8 ⩽ J e ⩽ 10 10 A m −2 , which overlaps quite well with the current density range explored in our work. Our results thus suggest that a small part (up to ∼ 15% here) of the current extracted at high current density regimes in experimental arrangements similar to Spindt (1992) could be assigned to a thermal coupling by the bulk between neighbouring emitters. In final words, it is interesting to note that, beyond thermal coupling, taking into account thermal effects to roughly estimate the pre-breakdown voltage significantly changes the studies aimed at optimizing spacing in FEA. Indeed, computing the temperature is a reasonable way to evaluate the maximum applied electric field under which an FEA electron source can operate without failure. In addition, it is important to take into account the gain in current due to the temperature. From these considerations, a densely packed distribution of sharp tips emitting in the thermo-field regime under a high applied electric field, that compensates for the screening but do not causes breakdown, can lead to a higher emitted current per unit surface than an electrostatics-only optimized distribution of the same emitters.
Taking an example from the present work, the emitted current per protrusion for f = 10 under 0.4 GV m −1 in the array arrangement with spacing d pb = 5.8 µm (table 3) is the same, I = 17 mA, than that of an isolated protrusion of f = 10 under 0.16 GV m −1 (figure 5). In terms of current per unit surface however, an array of 'isolated' protrusions with f = 10 requires at least a spacing of 50 µm (to have δ a ⩽ 1%) which gives I/d 2 = 6.8 × 10 6 A m −2 , while the array arrangement, being denser, delivers a higher current per unit surface, I/d 2 = 5.1 × 10 8 A m −2 , despite a strong apex electric field loss δ a = −60.4%.

Conclusion
This modelling work deals with thermo-field electron emission from idealized protrusions. It focuses on hemiellipsoids with aspect ratios f = H/R from 10 to 1 in order to model protrusions ranging from typical electron source emitters with notable apex FEFs, β a ∼ 50, to simple surface asperity with low apex FEF, β a ∼ 3. A size of H = 10 µm was chosen, but the results are consistent over variations of a few micrometers.
The transition over the applied electric field from pure field regime to thermo-field regime is characterized for the given set of protrusions. Both regimes are distinguished by a temperature contribution to the current respectively below and above 1%. At pre-breakdown voltage, the self-heating of protrusions demonstrates significant gains in current, ranging from +60% for f = 10 to +110% for f = 1. These maximum values outline the influence of temperature when dealing with breakdown problems in electron sources or high-voltage devices. It also highlights the need to account for thermal effect in the modelling studies aimed at optimizing spacing in FEA. In addition, it is observed that thermal contribution is higher for rounder protrusions, with smaller apex FEF. Through the analysis of local variable distributions, this result is attributed to the protrusions geometry: for smaller aspect ratios, an overall smaller local electric field spread over a wider emission area is responsible for a higher contribution of temperature to the current.
Colourmaps of the temperature distribution at equilibrium are shown at the pre-breakdown voltage for isolated protrusions. They give insights into the way protrusions of different shapes-here, different aspect ratios-heats up all over their volume until the melting point is reached. Besides, in line with previous works, our results show how a cooling Nottingham effect at the protrusion apex affects the equilibrium temperature distribution. This causes the highest temperature to be found not at the apex of the emitter but slightly below. The melting temperature being reached underneath the emitter surface, this could be related to an explosive breakdown phenomenon as reported in several experimental studies.
These temperature maps obtained at pre-breakdown also demonstrate significant heating of the cathode bulk around their bases. This observation supports the possibility of a thermal coupling between close protrusions via the bulk material on which they are located. In the thermo-field regime, modelling results for multiple protrusions in interaction demonstrate the effect of a thermal coupling from the base, up to a contribution to the current ε I of 13%. This maximum value is obtained for the smallest aspect ratio considered, f = 1. Indeed the heat around their bases was observed to spread wider for protrusions with lower aspect ratios. It should be noted, however, that the value of the pre-breakdown field for such dome-like protrusion ( f = 1) is out of the experimental range. For sharper protrusions, emitting at lower applied electric fields, the heat spreads over a smaller distance. The self-heating process for protrusions interacting at these smaller distances is drastically reduced by the electric field loss around the apex, therefore preventing any thermal coupling at their bases. The overall weak influence of thermal coupling on the emission is in fact related to the prevalence of the electrostatic screening when protrusions are densely distributed.
Nevertheless, for specific configurations achievable in FEA experiments, the electrostatic screening can be compensated by a higher applied electric field. A parametric study at a borderline applied electric field of 0.4 GV m −1 demonstrates noticeable thermal couplings (ε I ⩾ 5%) for protrusions sufficiently enhancing the electric field: f ⩾ 5, β a ⩾ 20. An upper value of ε I = 14.5% is observed for an array arrangement of protrusions with f = 10 and a spacing distance of 5.8 µm.

Appendix A. Emission current density calculation
The emission current density is obtained by integrating over all normal-energy ε n the product between the supply function N(ε n , T) and the transmission probability D(F, ε n ): The normal-energy writes ϵ n = ϵ C + mv 2 z /2 where v z is the velocity in the direction normal to the metal-vacuum interface and ε C is the energy at the bottom of the conduction band, with the zero-energy reference set at the vacuum level.
The supply function can be written as N(ϵ n , T) = 4πm h 3 k B T ln 1 + exp where ε F is the Fermi level. It is obtained from the Fermi-Dirac distribution in the framework of Sommerfeld description of metal.
The transmission is assumed to occur through a Schottky-Nordheim (SN) barrier whose potential E(z) along the normal coordinate z writes with its maximum value being E max = − e 3 F 4πε 0 (A4) and the transmission probability is obtained using the Kemble formalism (Kemble 1935) for a SN barrier. However, in accordance with the prescription of Murphy and Good (1956), the mathematics split into two cases at the normal-energy ϵ * n = E max / √ 2 (above the top of the potential barrier). For normalenergy above ϵ * n the transmission probability is set to 1, while for ϵ n ⩽ ϵ * n the formula writes In our simulation routine, all integrals, including the elliptic integrals (Press et al 1992) involved in v( y), are computed numerically.