Particle-in-cell simulations of high frequency capacitively coupled plasmas including spatially localised inductive-like heating

High frequency (HF) capacitively coupled plasmas (CCPs) are ubiquitous, having several industrial applications, especially in the semiconductor industry. Inductive heating effects within these plasmas play an important role and therefore understanding them is key to improve industrial applications. For this purpose kinetic research, using particle-in-cell (PIC) codes, offers significant opportunity to study, and improve, industrial plasma processes that operate at the atomic level. However, PIC codes commonly used for CCPs are electrostatic and thus cannot be used to simulate electromagnetically induced currents. Therefore we have developed EPOCH-LTP, a 1D PIC code with a current heating model, that enables the simulation of inductive heating effects in HF CCPs. First simulation results, from an HF CCP (60 MHz) operated at 1 mTorr of argon, show that inductive currents couple most of their power to the electrons at the interface between the bulk plasma and the sheath. Furthermore, the simulation of a dual-frequency CCP, where a HF inductive current and a low-frequency (LF) voltage waveform at 400 kHz are applied, have shown a synergy between the HF and LF waveforms that increase the inductive heating rate.


Introduction
Low pressure capacitively coupled plasmas (CCPs) play an important role in industrial and biomedical applications [1,2].In the semiconductor industry CCPs are used for atomic layer etching and deposition processes that are essential for the production of integrated circuits [3,4] where the ions in the plasma are directed towards a wafer.At present, the precision has reached the atomic level [5,6] and so there is a demand for a high degree of control of the ion flux, ion energy and ion angle of incidence.
The application of relatively high-frequency (HF) driving voltage in the megahertz range, ≳60 MHz, are of significant interest to increase the plasma generation efficiency [7] and, when combined with a low frequency (LF) voltage waveform, allows separate control of electrons and ions [8,9].With higher frequency waveforms the sheath impedance decreases and thus enables a larger degree of electron power coupling in exchange for lowering the ion energy and electrode voltage [10][11][12].However, with increasing driving frequencies, electromagnetic effects become significant if the excitation wavelength λ and the plasma skin depth δ are of the order of magnitude of the plasma source [13,14].In such cases there is a transition from capacitive to inductive heating [13,14] that causes transverse inductive power coupling [15][16][17][18][19] and strong plasma inhomogeneities [20][21][22][23], which can have a significant impact on the ion fluxes [24,25].While there has been progress in studying impact of electromagnetically induced currents on the ion kinetics [10,26], further research is needed to understand these phenomena in detail.
It is therefore of interest to develop a detailed understanding of how inductive heating can impact the operation of HF CCPs, and in particular with respect to the distribution of ion fluxes, energy and angle of incidence at the electrodes.For this purpose, we have developed EPOCH-LTP, a 1D3V electrostatic particle-in-cell (PIC) code with an inductive heating model.EPOCH-LTP is the low temperature plasma spin-off of the highly parallelised code EPOCH [27].The inductive heating model is based on the work of Meige et al [28] that has been used before in PIC simulations for inductively coupled plasmas [28][29][30][31].The inductive heating model is a heating mechanism that imposes a perpendicular oscillating current that replaces the inductive currents present in HF CCP.We have developed Meige's model further by improving its spatial resolution and have implemented it in EPOCH-LTP.This implies a fundamental change to the already existing inductive heating model, and thus the goal of this paper is to understand the behaviour of this new model and how it could be used to investigate inductive heating effects in HF CCPs.
For this investigation we propose a HF (60 MHz) CCP operated in argon at low pressure (1 mTorr = 0.13 Pa) as a case of study.First, the methodology is described in section 2, and contains a description of the recently developed PIC code EPOCH-LTP in 2.1, the collision model in 2.2 and code validation in 2.3.The inductive heating model is described in 2.4, and a detailed description of the HF CCP simulation setup is found in 2.5.The simulation results are described in section 3, where the HF CCP case of study is first characterised in 3.1, followed by a parametric investigation of the inductive current amplitude in 3.2, and an applied direct current (DC) bias in 3.3.These studies build sufficient confidence to present, in 3.4, the simulation results for a dual-frequency driven CCP with an HF inductive current and an LF voltage waveform.Finally, a summary, conclusions and future work is proposed in section 4.

Method
The use of PIC models for kinetic investigations of CCPs is well established within the research community [32][33][34].PIC models have been used to investigate HF CCPs [7,10,35,36] however, they implement an electrostatic field solver [34,37] and therefore the investigation of electromagnetic effects is restricted.To overcome this limitation Meige et al [28] implemented an inductive heating method in a 1D PIC model that has enabled kinetic investigations of inductively coupled plasmas [29][30][31].Based on this, we have implemented a similar inductive model into EPOCH-LTP, a 1D3V electrostatic PIC code, that enables kinetic simulations of inductive heating effects in HF CCPs.This section first describes EPOCH-LTP, followed by a detailed description of the numerical implementation of inductive currents, and finally describes the simulation configuration for HF CCPs.

PIC code EPOCH-LTP
In this work we make use of the recently developed PIC code EPOCH-LTP, a 1D3V highly parallelized code developed from the parent code EPOCH, a PIC model for high energy density plasmas [27].EPOCH-LTP solves Poisson's equation for the electric potential with a finite-difference method [37], integrates particles in phase-space using a Boris scheme [34], particle-grid interpolation using triangular (for particles) and squared (for grid cells) shape functions [34], and models electrodes as perfect absorbing walls.For operation in pure argon, we consider electrons, ions (Ar + ) and neutrals where the latter are treated as an spatially uniform background field with constant density and temperature.

Collision model and reaction chemistry
The collisions between neutral and charged particles are modelled implementing a null-collision Monte Carlo collision (MCC) method [38,39].The relevant collision processes are shown in table 1.They are elastic scattering (reaction number 1), excitation (2) and ionisation (3) [40,41] for electronneutral collisions and ion-neutral elastic scattering (4) and charge-exchange (5) [42].All collisions are assumed to be isotropic, excited states are not tracked, i.e. excitation processes act only as an electron energy loss mechanism, and the residual energy resulting from electron impact ionisation is equally divided between the two outgoing electrons.

EPOCH-LTP validation
EPOCH-LTP has been tested against established PIC simulation results [35,[43][44][45] and experimental work [46] to ensure  The electron and ions density profiles, shown in figure 2(a), show a symmetric distribution with a quasi-neutral central region, where n e ≃ n Ar + ∼ 7 × 10 15 m −3 .The bulk plasma is surrounded by sheaths where a faster density drop for electrons than for ions is observed, which is consistent with expectation.The plasma potential, in figure 2(b), also presents a symmetric distribution with a flat potential around ∼112 V where the plasma is quasi-neutral and drops across the sheath to match the 0 V time averaged potential at the electrodes.The spatially resolved profiles of plasma potential and density both show good qualitative agreement with Donko et al [44] and with similar simulation problems in [1,33,45,47].Moreover the electron energy probability function (EEPF) and the ion flux energy distribution function (IF-EDF), shown in figures 2(c) and (d) respectively, also show good agreement.The EEPF shows expected results for low pressure Ar CCP [44,48] and the IF-EDF shows the characteristics peaks caused by ionneutral charge-exchange collisions as the ions cross the sheath [49], demonstrating good agreement with experimental results [46].
Two further simulation studies are undertaken, which suggest that EPOCH-LTP is operating consistently with expectation, are briefly described here.Further information on the corresponding results is provided in the electronic supplementary information.First, we study the benchmark problems from Turner et al [43] for CCP operating in helium.These are successful and bring EPOCH-LTP in line with other existing and well-established PIC models.The spatial profile of the ion density of these benchmark problems can be assessed quantitatively with their corresponding χ 2 distribution functions provided in [43].The χ 2 values obtained for problems 1 to 4 with EPOCH-LTP are within the 95% confidence interval, which suggests that EPOCH-LTP is operating correctly.Please refer to figure S1 of the electronic supplementary information for more details.In a second study, we consider the argon CCP investigation from Schulenberg et al [46], and for a pressure variation (1-10 Pa) observe good agreement with their experimental and simulation results.These tests also show good agreement in the electron peak density and in the IF-EDF, as shown in figures S2 and S3 of the electronic supplementary information.

Inductive heating method
The inductive heating caused by electromagnetically induced currents is modelled applying an external current perpendicular to the simulation domain.Note that because this heating mechanism is externally imposed and is applied in the context of an electrostatic field solver, and thus electromagnetic effects, such as skin depth, are not expected to be observed.
This model was first developed by Meige et al [28] and has previously been used in the context of inductively coupled plasmas [29,31] and space propulsion [30] research.The model is based on Ampere-Maxwell's law, which describes that the total induced current density is composed of a displacement current caused by time variations in the electric field E y (x, t), and by conduction current caused by transport of electric charges.Note that ε 0 is the vacuum permittivity, E y the electric field perpendicular to the simulation domain, and the subscript s refers to a given particle species that carries an electric charge q s , has a perpendicular particle flux Γ y,s = n s vy,s , number density n s and a mean velocity vy,s in the perpendicular direction to the axial distance between the electrodes.The inductive heating model consists of imposing a current density, J(x, t), over an arbitrary spatial extent and computes the resulting perpendicular electric field, E y , generated by the displacement current J disp .The induced E y field is then used in the integration of the particle equations of motion.
Regarding the time resolution, the numerical implementation is executed for each simulation cycle, t → t + ∆t, and the Courant-Friedrichs-Lewy (CFL) time restriction is considered to ensure sufficient accuracy.
Regarding the spatial resolution, unlike in Meige et al [28] where E y is calculated globally over the whole inductive source, the method described here is executed individually on each grid cell i, of width ∆x.This improves spatial resolution, and enables us to account for interactions between the plasma and inductive currents that can vary with respect to spatial location between the electrodes.
The numerical implementation, for each time step t → t + ∆t and for each grid cell i, is as follows.First the plasma conduction current is computed by making use of (3) where W s is the super-particle weight factor.Secondly, the perpendicular electric field time differential is computed using equations ( 1) and (2) Finally, the electric field is integrated in applying an implicit Euler method The time integration of E y requires an initial condition at t = 0 which, assuming no initial induced field, is set to The inductive heating method has been tested reproducing four simulation cases, labelled A to D, from Meige et al [28].These test problems consist of a 1D inter-electrode gap of 10 cm, containing an argon at 1 mTorr (0.13 Pa) and 297 K, and grounded electrodes, as shown in figure 3. The inductive current source, Figure 3. Simulation setup for an argon plasma driven with a HF inductive current, J(x, t).J(x, t) is applied in the y-direction between x min and xmax (green shaded area). is driven at f = 10 MHz and J 0 = 100 Am −2 and placed between x min and x max in space.The simulations run for 1000 cycles, with a time step resolving the CFL condition for 100 eV electrons and a grid width that ensures the Debye length is resolved, ∆x < λ D .The reaction set included is as described in section 2.2 and the simulation parameters are listed in table 3.
Figure 4 shows the resulting spatial profiles for Ar + and the plasma potential, where the spatial extents for the application of inductive heating corresponding to case A-D are shown in figure 4(a).The resulting plasma density and potential for these tests are shown in figures 4(b) and (c), respectively.The inductive heating method is performing as expected, generating an homogeneous plasma independent from the location and extent of the inductive source.Even though the heating Steady state parameters for an argon plasmas test problem described in [28].In (a) the spatial extent of the inductive source J(x, t), and the corresponding spatial variations of (b) ion density and (c) plasma potential.Simulation parameters listed in table 3. source is not symmetrically located, the resulting spatial profile is quasi-symmetric.This is because electron mean-freepath is larger than the reactor characteristic length, λ mfp,e > L. and thus ionisation processes are non-local and are found in the entire simulation domain.
The plasma potential, ϕ p ∼ 17 V, is as expected for a plasma generated between floating walls with collisionless sheaths [1].The plasma potential profiles are in good agreement with the results in [28] that are expected to be near the ionisation energy of argon, E ion .There is however a difference between applying the inductive source at the edge (cases A and B) or in the bulk plasma (cases C and D), that lead to different plasma densities.which results in different peak values of the plasma density: ∼2 × 10 16 m −3 and ∼3 × 10 15 m −3 , respectively.
On the one hand, when the inductive source is placed at the bulk plasma (C and D) the density results, shown in figure 4(b), are in good agreement with results in [28].Note that the peak density in case C is slightly higher than in case D because longer inductive sources couple more power to electrons.
On the other hand, cases A and B show densities that are an order of magnitude higher than in cases C and D. This is caused by enhanced power absorption within the sheath that underpins ionization and thereby the plasma density.In these cases the penetration of the heating source into the bulk plasma does not play an important role and therefore the density for cases A and B are similar.This behaviour, discussed in detail in section 3.1, is not observed in [28] and this can reasonably be attributed to the relatively high spatial resolution that is implemented for the purposes of this work.
Previous work using the inductive heating method did not observe enhanced power deposition within the sheath because the resulting displacement currents were averaged over the entire, instead of being discretized into grid cells, and therefore regions with high displacement currents were smoothed out [28].Nevertheless, these preliminary simulations show that the inductive heating method implemented in EPOCH-LTP is consistent with previous literature [28][29][30][31].

Simulation setup for HF CCPs
CCPs are typically generated in cylindrical reactors, with circular electrodes of radius R separated a distance L in the axial direction.Assuming that the plasma is axially symmetric, the model simulates the axial direction at a given radial position, r, with the inductive current pointing in the radial direction.
In HF CCPs standing wave effects (SWEs) are responsible for the transition from capacitive to inductive heating in the radial direction [13,14].For weak SWEs (intermediate R values), capacitive heating is maximum at the centre (r → 0) and decays in the radial direction, and inductive heating is near zero at the centre and peaks at outer radial positions (r → R) [14].This cause non-homogeneous power deposition, however at low pressure, ⩾ 30 mTorr, the electron mean free path, λ mfp,e , is typically larger than the reactor size, λ mpe,e ≫ R, and so the density profile is spatially homogenous as it is determined by diffusion rather than by power deposition [14].
Besides, the sheath impedance in CCPs is inversely proportional to the driving frequency [1, ch 11] and, at a certain threshold frequency, the impedance reaches a minimum [10,11].At the limit of minimum sheath impedance, the radio frequency potential amplitude can be neglected and the plasma reactor walls can be defined as floating walls [10].
Based on these considerations, it is reasonable to model outer radial positions of HF CCPs operated at low pressure that are sustained only by an inductive transverse current.Therefore, the simulation setup consists of two grounded electrodes containing argon at 1 mTorr (0.13 Pa) with a perpendicular inductive heating source, as sketched in figure 3. Note that the simulation setup is designed to study the effects of inductive heating in isolation and therefore other important effects such as secondary electron emission.
The heating source is defined as with f = 60 MHz, and amplitude that is applied uniformly from the left electrode until about x S = 2 cm, where a hyperbolic tangent function avoids a sharp discontinuity, as shown in figure 5. Further simulation parameters are listed in table 4, and the collisions included are as described in section 2.1.

5.
Spatial profile of the inductive current amplitude J 0 (x).The grey dotted line denotes the position of x S , where J 0 falls to half its amplitude.The simulation domain is split into 2500 cells, ensuring enough resolution for resolving the plasma Debye length, i.e. λ D ≫ ∆x.The simulations run for 7000 cycles with a time step that satisfies the CFL condition for electrons at 100 eV.The simulations converge within 3000 cycles and steady state and phase resolved data is extracted over the last 500-1000 cycles.
The simulation parameters described in this section are used in the following sections unless otherwise explicitly stated.

Characterization of inductive heating effects in HF CCPs
Kinetic effects of inductive heating in HF CCPs have been characterised with the results of the simulation setup described in section 2.5.Steady state results are presented in figures 6 The steady state results, in figure 6(a), show an homogeneous and symmetrical plasma density distribution, with n e ∼ 10 16 m −3 at the midplane.As observed in section 2.4, the profile is symmetric because of the non-local nature of electron impact ionization processes with respect to the heating source, λ mfp,e > L, and the plasma potential, in figure 6(b), is determined by the floating wall boundaries, ϕ p ∼ 17 V.The EEPF, infigure 6(c), shows the electron population crossing the midplane (x = 5 cm) from left to right, i.e. electrons moving away from the inductive source.Most of the electrons are found below the ionisation threshold, ε < ε ion .These electrons are in thermal equilibrium following a Maxwellian distribution with T e = 3.04 eV.The electrons that are at higher energies, ε > ϕ P , present populations several orders of magnitude lower.This indicates that electrons in the bulk plasma either rapidly dissipate energy via electron-neutral collisions or are lost across the sheaths.The IF-EDF at the left electrode (x = 0 cm), figure 6(d), shows a single peak at ε ∼ ϕ p , indicating that the ion acceleration across the sheath is collisionless and is not significantly affected by displacement currents, which is in accordance with the experimental observations of [25].
The inductive power, P e,y = E y J e,y , is mostly coupled to the electrons due to its lower inertia.This is particularly high in the sheath, as shown in figure 7, and thus the inductive power absorption in this region determines the plasma formation.The power coupling in the bulk plasma, compared to the sheath, is negligible and therefore the penetration of the inductive current into the plasma does not play an important role in the plasma formation.
More detail about the high power coupling at the sheath can be obtained from the phase-resolved data near the left electrode (0 ⩽ x ⩽ 5 mm), as shown in figure 8.The electron density, in figure 8(b), shows that the plasma sheath is steady in time and thus not perturbed by the inductive current, displayed in figure 8(a).The sheath is ∼1 mm wide and, as expected, presents a sharp density drop.It is this density drop that causes the high inductive power coupling in the sheath.
Spatially localised values of the plasma density play an important role in power absorption because the response of the plasma to the inductive current J(x, t) depends strongly on the number of particles present per unit volume.The interaction between the inductive current J(x, t) and the plasma can be described as the plasma conduction current J cond constantly adjusting to the inductive current.This adjustment is carried out by the displacement current J disp which generates the electric field E y needed to redirect the particle velocity.In case there are many particles present per unit volume, i.e. high density, the velocity redirection required at each particle is small and therefore a smaller E y is generated.However, when there are fewer particles (lower density), the change in particle velocity is more important and higher E y is required.It is in this case when the induced E y field is high and the inductive power absorption P e,y becomes maximum.When the limit x → 0 mm is approached, E y continues to grow, but n e is zero and so is J e,y , and therefore the absorbed power decays to zero.The density in the bulk plasma, x > 1 mm in figure 8(b), is high enough for J e,y to respond rapidly to variations in J(x, t).Consequently small E y fields are induced.The temporal variation of J e,y , shown in figure 8(c), presents an oscillating behaviour very similar to J(x, t), and E y , in figure 8(d), is small.This is confirmed in figure 9 where J e,y (red dashed line) and E y (blue dashed line) at x = 4 mm are plotted in time.These curves show that J e,y in the bulk plasma responds quickly to J(x, t) with little aid from E y .
The density in the sheath (x < 1 mm) is much lower than in the bulk plasma and thus variations in J(x, t) have a bigger impact on the particle velocity.Figure 8(d) shows that stronger E y are formed and accelerate the electrons until J e,y matches J(x, t).The required velocity change is large enough that the current match is not instantaneous, but a dephase is observed between the inductive and the conduction current.The dephase is inversely proportional to the density, as can be seen in 8(c) where J e,y reaches its maximum later in time the closer it is to the electrode, where n e → 0. This trend also correlates with the induced E y , which is stronger with decreasing density, as shown in figure 8(d).Curves for J e,y and E y at x = 0.5 mm are plotted in figure 9, showing that J e,y not only presents a delay with respect to J(x, t) but also a slight amplitude overshoot that is caused by the electron inertia.It is also observed that E y is much larger than in the bulk plasma and therefore the inductive power coupling is much more efficient in the sheath, as shown in figure 8(e).These results are in agreement with results presented in [22,23], where large radial electric fields and power absorption rates near the sheath were also observed.

Current amplitude variation
This section investigates how variations in the inductive current amplitude, J 0 in equation ( 9), influence the plasma.A set of simulations have been undertaken where J 0 has been varied from 1 to 500 A/m 2 .The super-particle weight has also been adjusted between W = 2 × 10 7 −10 10 in order to find a balance between ensuring good statistical values [34,50] and a reasonable computation performance.The remaining simulation parameters are as described in table 4 and results are shown in figure 10.Variations of J 0 are reflected in the peak plasma density and the peak electron power absorption, as shown in figure 10.The remaining plasma parameters are quantitatively similar to those shown in figures figures 6 and 7, i.e. uniform spatial distribution, plasma potential ϕ p ∼ 17 V, and EEPF and IF-EDF with similar features.
The electron power absorption is dominated by inductive heating at the sheath, as observed in section 3.1.Increasing J 0 however leads to larger induced displacement currents, and thus to larger induced E y fields.Therefore higher electron power coupling rates are found in the sheath, as shown in figure 10(a), allowing for higher ionisation rates and a correspondingly denser plasma, as shown in figure 10(b).
Therefore, J 0 is a parameter that controls the power coupled to the plasma, with a direct impact on the electron power absorption rate and the plasma density.

DC bias variation
This section investigates the effects of a DC voltage on a plasma sustained by HF inductive currents.For this, the left electrode is powered with a DC voltage, as sketched in figure 11, and has been varied between V DC = 0 V and 300 V.The resulting plasmas present an asymmetry with a wider sheath at the powered electrode, as shown in figure The plasma potential is ϕ p ∼ 17 V, as in previous cases, however the sheath width grows with V DC because larger potential drops must be shielded at the powered electrode, as shown in figure 12(b).This has an impact on the maximum energy that electrons can gain, and is reflected in the high energy tails of the EEPFs, shown in figure 12(c).Because of the wider sheath electrons are able to penetrate deeper into the sheath, where the displacement currents are larger, and thus gain more energy from the inductive heating.This results in a population increase for ε > ϕ P although it does not significantly affect the plasma production rate as these electrons are more likely to escape the plasma through the opposite sheath before undergoing any collision.The IF-EDFs, shown in figure 12(d), present peaks at ε ∼ ϕ p + V DC and with lower peak population for increasing V DC .This decrease in the peak population is caused by the broadening of the peak, i.e. ions crossing the sheath gain or lose a small fraction of energy because of the large displacement currents present near the electrode.The higher the V DC applied the longer the ion exposure to high inductive fields while crossing the sheath and thus the wider the IF-EDF peak broadening is.However, the total ion flux at the electrode is equal for the four cases shown in figure 12(d) and is therefore independent from V DC as expected.
The electron power absorption P e,y , as shown in figure 13, presents a peak in the sheath as observed before.The peak value is constant, as expected because this depends on J 0 .However its location varies with V DC , shifting inwards with increasing V DC .This shows that the most effective power absorption location is found at the interface between the sheath and the bulk plasma, where the density is relatively low, but non-zero.In summary, the application of a DC bias introduces an asymmetry to the plasma, with a wider sheath at the powered electrode that shifts the electron power absorption peak.

inductive heating with low-frequency oscillating potential
In this section we investigate how a LF voltage waveform interacts with a HF inductive current.The simulation, as sketched in figure 14, is configured to apply an oscillating voltage source to the left electrode at a frequency f L = 400 kHz and amplitude V 0 = 300 V. Note that the −1 term introduces a negative DC bias V DC = −V 0 .
The remaining simulation parameters are as described in table 4. The simulation results are presented in figures 15  and 16 and show steady state and phase resolved plasma parameters, respectively.The resulting plasma presents an asymmetry, as shown in figure 15(a) that is caused by the DC bias term in V s (t).The plasma potential, as shown in figure 15(b), is ϕ p ∼ 17 V, as in the cases presented in previous sections.The EEPF, shown in figure 15(c), is very similar to the distributions shown for cases with a DC bias applied.This, as described in section 3.3, is caused by the presence of a wider sheath.The IF-EDF presents a two-peak distribution as typically found in low frequency driven CCPs [1, ch 11].
The phase-resolved results show that the LF voltage waveform, black line in figure 16(a), is responsible for the expansion and collapse of the sheath.The electron density, presented in figure 16(b), shows that the sheath is completely expanded when V s (t) is at its minimum (t = 0 and 2.5 µs), and collapses when V s (t) = 0 V (t = 1.25 µs).The inductive power absorption peaks at the sheath edge, as described in section 3.3, and therefore a high power absorption trace that follows the sheath motion is observed in figure 16(c).
is especially high when the sheath is collapsed, as shown in figure 16(a) (red solid line) where P T e,y is integrated between the limits x = 0 mm and x = 8 mm.When the sheath shrinks (0 < t < 1.25 µs) it allows a plasma expansion that lowers the plasma density near the sheath.This widens the extent over which the displacement current is important, and thus inductive power absorption becomes more effective.When the sheath is fully collapsed (t ∼ 1.25 µs), the extension over which the power absorption is important is at its maximum and, it is in this time window when electrons gain most of their energy.The mean electron energy, shown in figure 16(d), shows that electrons gain more energy while the sheath is collapsing and peak when the sheath is fully collapsed.This causes an increase in the ionisation rate, by a factor of about ∼5, as shown in figure 16(e).When the sheath begins to expand (t > 1.25 µs), the opposite effect occurs, i.e. the plasma contracts, the density increases and the effective power absorption length narrows and thus the inductive current couples less energy.
Therefore, the interaction of a HF inductive current and a LF voltage demonstrates interesting behaviour as dynamics at very different time scales are able to create synergies that enhance the inductive power coupling and plasma generation.This is clear when comparing results with the V DC = −300 V case in section 3.3.The presence of a LF voltage increases the total power absorption by 30% and the plasma density increases by approximately 20%.
These results can serve as a basis for further investigation of kinetic inductive heating effects using 1D electrostatic PIC codes.Apart from the parameters studied in this work, there is a wide range of parameters that would be useful to investigate to understand inductive heating in more detail, such as the driving frequency, the gas pressure or the spatial distribution of the inductive source.Finally, it would be worth to study the interaction with LF waveforms and how synergies could be optimised.

Conclusions
In this paper we present EPOCH-LTP, a high parallelized 1D/3 V PIC code including an inductive heating model.This model is proposed to study kinetic inductive heating effects in HF CCPs.EPOCH-LTP has been used to simulate a CCP operated in argon at 1 mTorr and driven at 60 MHz.The simulation results show that transverse inductive currents in the sheath dominate the inductive heating and thus the plasma formation.The plasma density is directly proportional to the size of the inductive current.The application of a DC bias shifts the location of the inductive heating peak.When the HF inductive current is combined with a LF voltage waveform the inductive heating is controlled by the sheath dynamics and this improves the plasma generation efficiency.This is an interesting result as dynamics at very different time scales are able to create synergies that enhance the inductive power coupling and thus the plasma density.

Figure 4 .
Figure 4.Steady state parameters for an argon plasmas test problem described in[28].In (a) the spatial extent of the inductive source J(x, t), and the corresponding spatial variations of (b) ion density and (c) plasma potential.Simulation parameters listed in table 3.

Figure 6 .
Figure 6.Spatial variation in the (a) electron and ion density and (b) plasma potential for the simulation setup described in section 2.5, parameters listed in table 4, and sketched in figure 3. Corresponding (c) EEPF at midplane (x = 5 cm) and (d) IF-EDF at the left electrode (x = 0 cm).

Figure 7 .
Figure 7. Spatial variation in the steady state electron power absorption for the simulation setup described in section 2.5.

Figure 8 .
Figure 8. Phase resolved (a) normalised inductive current J(x, t), (b) electron density ne, (c) perpendicular electron current density Je,y, (d) perpendicular electron field Ey, and (e) perpendicular power absorption Pe,y for the simulation setup described in section 3.1.The vertical axis in panels (b)-(e) present the first 5 mm next to the left electrode (x = 0 cm), as sketched in figure 3.

Figure 9 .
Figure 9. Phase resolved inductive current density (black), electron current density (red), and induced electric field (blue) at different locations.Current densities and electric fields are normalised with J 0 = 100 A m −2 and E 0 = 8 × 10 3 V m −1 , respectively.

Figure 10 .
Figure 10.(a) Peak electron power absorption in the sheath and (b) electron density at the simulation midplane (x = 5 cm) under variations of J 0 for the operation conditions listed in table 4.

Figure 11 .
Figure 11.Simulation setup for an argon plasma driven with a HF inductive current J(x, t) and a powered electrode with a direct current (DC) voltage V DC .

Figure 12 .
Figure 12.Spatial variation in the (a) electron density and (b) plasma potential for a plasma driven by HF inductive current with a DC voltage, as described in section 3.3 and sketched in figure 11.Corresponding (c) EEPF at the midplane (x = 5 cm) and (d) IF-EDF at the powered electrode.

Figure 13 .
Figure13.Spatial variation in the electron power absorption under different applied V DC bias voltages.The colour code is as in figure12.

Figure 14 .
Figure 14.Simulation configuration for an argon plasma driven with a HF inductive current J(x, t) and a LF voltage Vs(t).

Figure 15 .
Figure 15.Spatial variation in the (a) electron and ion density and (b) plasma potential for an argon plasma driven by a HF MHz) inductive current and a LF (400 kHz) voltage waveform.Operation conditions are listed in table 4 and the simulation configuration is sketched in figure 14.Corresponding (c) EEPF at midplane (x = 5 cm) and (d) IF-EDF at the powered electrode (x = 0 cm).

Figure 16 .
Figure 16.Phase resolved (a) LF voltage waveform Vs(t) (black) and total electron power absorption P T e,y (red), (b) electron density ne, (c) electron perpendicular power absorption Pe,y, (d) electron mean energy εe, and (e) ionisation rates R i for a plasma driven by a HF inductive current and a LF voltage waveform, as described in section 3.4 and sketched in figure 14.The vertical axis in the heatmap plots present the first 8 mm next to the powered electrode.

Table 1 .
Collisions implemented in Monte Carlo collision (MCC) model.εexc and ε ion are the excitation and ionisation energy thresholds, respectively.

Table 3 .
[28]lation parameters for an inductively heated argon plasmas proposed in[28]and sketched in figure3.The corresponding results are shown in figure4.

Table 4 .
Simulation setup for a HF CCP.