TurnerT Case 2 (TurnerT.pre)


capacitively coupled plasma, CCP, discharge, steady state, TurnerT

Problem Description

In this example we demonstrate VSim’s ability to simulate capacitively coupled plasmas, using the benchmark cases of Turner et al. [TDD+13]. Turner’s work documents the successful benchmarking of five independently developed particle-in-cell codes (not including VSim) for four different capacitively coupled discharges at various background pressures.

Here, we consider the second of the Turner scenarios, though the input file can be readily modified to simulate the others. In addition to being able to accurately reproduce the Turner results, VSim can also employ physics-based initialization methods to enable more rapid convergence of the simulations to their steady-state. The use of such methods will also be explained below.

This simulation can be performed with a VSimPD license.

Opening the Simulation

The TurnerT example is accessed from within VSimComposer by the following actions:

  • Select the NewFrom Example… menu item from the File menu.
  • In the resulting Examples window expand the VSim for Plasma Discharges option.
  • Expand the Capacitively Coupled Plasmas (text-based setup) option.
  • Select Turner case 2 (text-based setup) and press the Choose button.
  • In the resulting dialog, create a New Folder if desired, and press the Save button to create a copy of this example.

The basic variables of this problem should now be alterable via the text boxes in the left pane of the Setup Window, as shown in Fig. 476.

image 1

Fig. 476 Setup Window for the TurnerT example.

Input File Features

The basic physics of this simulation is a balance between collisional processes and wall losses; a one-dimensional box of length 6.7 cm contains neutral helium gas at room temperature (300 K) and density 3.21e21 1/m^3. The gas is weakly ionized, resulting in a population of free electrons and singly ionized helium atoms at density 5.12e14 1/m^3. The helium ions are also at room temperature, while the electrons are considerably hotter (30,000 K). The left wall of the box is grounded, while the right wall oscillates with a bias voltage of 200 V at frequency 13.56 MHz.

Charged particles are lost upon collision with the wall, and are replenished by ionization of the background neutral gas by the hot electrons; the latter process repopulates both the electrons and helium ions in the plasma (the background neutral gas is treated as an infinite source). Plasma sheaths form near the walls, containing electric fields which are strong relative to those elsewhere in the plasma; the particle density profiles adjust in response to the fields in the sheath. The sheath transit time, for ions, is much longer than the period of the oscillating potential; thus, multiple RF cycles occur while an ion crosses the sheath. A steady state is attained when the loss rate of particles to the wall comes into balance with the ionization rate for a particular profile shape.

In our initial run we are not going to model the full evolution of the discharge to its steady-state parameters; rather, we will explore the basic physics of the discharge and modify the simulation accordingly (with the aim of ultimately hastening convergence to this steady state, while exploring VSim capabilities).

Running the Simulation

To run the simulation, continue as follows:

  • Proceed to the Run Window by pressing the Run button in the left column of buttons.
  • The Number of Time Steps is 10240 (1/400 of the length of Turner’s run) for this example, and can be modified if desired.
  • Consider checking the ‘Run In Parallel’ box if desired (set a ‘Number of Cores’ appropriate for your VSim license).
  • To run the file, click on the Run button in the upper left corner of the right panel. You will see the output of the run in the right pane. The run has completed when you see the output, “Engine completed successfully.” A snapshot of the simulation run during execution is shown in Fig. 477.


10240 steps takes roughly ten minutes on a MacBook Pro in serial, and less than four minutes on a four-core Windows machine.

image 1

Fig. 477 The Run Window during execution.

Analyzing the Results

We are going to run a postprocessing script, computePtclNumDensity.py, which builds density profiles from the particle data generated by VSim, so that we can look at these profiles and their evolution. To do so, we do the following:

  • Click the Analyze Window
  • Check the Show All Analyzers box
  • Select computePtclNumDensity.py from the list of analyzers and click the Open button at the bottom of the left panel
  • Fill in the text boxes
    • The simulationName should already be filled in, but if it is not, type in the name of the .pre file without the .pre extension, which is TurnerT.
    • For the speciesName, type in ‘electrons’ without the quotes.
  • Click the Analyze button in the top right of the right panel; this will generate the electron density profiles.
  • Now replace ‘electrons’ in the speciesName box with ‘He1’, for the helium ions.
  • Click the Analyze button again to generate the ion density profiles.

The resulting data will be visualizable as electronsDensity and He1Density under the 1-D Fields Data View in the Visualize Tab.

image 1

Fig. 478 Analysis window for the TurnerT example.

Visualizing the Results

Now that we’ve got all of our data, let’s look at it.

  • Click the Visualize Window

After a brief moment the visualization options for this data should appear.

We’ll first look at the time evolution of some fundamental one-dimensional quantities. From the Data View pulldown menu on the top left, select History. The default view here should contain four plots, namely, the electron and ion currents to the left wall and the number of electrons and ions in the simulation. A number of notable physics effects can be seen here:

  • After a sharp initial decrease in the electron population, both ion and electron populations decline at approximately the same rate. This is not as apparent from the separate numElec and numIons plots, but clicking on the “Location” dropdown window in Graph 3 and selecting “Window 4” as the new rendering destination, places both ion and electron populations in the same plot. (Select “<None>” in the plot variable (the topmost menu) for both Graph 1 and Graph 2 to resize the electron/ion plot.) The initial decrease in electron population arises when rapid electron wall losses create a charge imbalance in the plasma and establish plasma sheaths near the walls. Thereafter, this charge imbalance is preserved and the transport of both electrons and ions to the wall becomes ambipolar. A history of the particle popluations can be seen in Fig. 479.
image 1

Fig. 479 The electron and ion populations versus time.

  • The electron wall currents are quasi-periodic. The oscillating potential drives the highly mobile electrons alternately into the left and right walls. In the plot variable menu, change “numElec” to “leftElecCurr” in Graph 3 and “numIons” to “rightElecCurr” in Graph 4. The impacts of the electron cloud on the left and right walls, and their phasing in time, can be seen in response to the potential oscillations. A history of the electron currents can be seen in Fig. 480.
image 1

Fig. 480 Electron currents on the left (green) and right (black) walls versus time.

  • The ion currents are non-periodic. Ions, being much heavier than the electrons, exhibit relatively little response to the oscillating potentials. In the plot variable menu, change the Graph 1 quantity “None” to “leftIonCurrent” and the location to “Window 4”, then change the Graph 2 quantity “None” to “rightIonCurrent” and the location again to “Window 4”. The ion currents do not have the quasi-periodic structure of the electron currents; rather, ions diffuse outward to the walls in response to the DC sheath potentials, which are established by the initial departure of electrons and may also be rectified by the RF. A history of all the particle currents can be seen in Fig. 481.
image 1

Fig. 481 Electron and ion currents on the left and right walls versus time.

  • Ion losses are negligible before the initial establishment of the sheath. Change the plot quantity in Graph 3 from “leftElecCurr” to “None”. Change the plot quantity in Graph 4 from “rightElecCurr” to “numElec” and the “Location” to “Window 3”. It is clear that the dominant loss of ions to the wall only begins after the initial decrease in electron population (which corresponds to the establishment of the sheath). A history of the electron population against the ion currents can be seen in Fig. 482.
image 1

Fig. 482 History plots showing the majority of the ion current to the walls only begins after initial decrease in electron population.

We can also look at the plasma sheath and the ensuing changes in density profiles directly. In the “Data View” menu at the top left of the Composer window, select “1-D fields”. The plot controls here are similar to those of the history window. Select “YeeStaticElecFldTrilinos_x” for the plot variable in Graph 1. Select “YeeStaticElecFldTrilinosPotential” for the plot variable in Graph 2. Select “electronsDensity” for the plot variable in Graph 3. Select “He1Density” for the plot variable in Graph 4, and select “Window 3” for the location of this plot. The evolution of the discharge in time can be viewed by moving the time slider below the plots.

A number of additional physics features can be seen:

  • Sheath effects are present. Regions of sharp potential variation, corresponding to strong electric fields, arise near the walls, but such fields are screened out in the bulk plasma. Moving the time slider, it is clear that this sheath behavior persists regardless of the phase of the oscillating wall potential.
  • Electron profiles are altered much faster than ion profiles. Both ions and electron profiles are initially constant (5.12e14 1/m^3), but by the time the first nontrivial dump file is produced (at time dumpPeriodicity * dt, approximately 1/3 of the way through the period of the first wall oscillation), electron-poor regions corresponding to the sheaths have already been established in the electron profile, while the ions have barely begun to respond to the presence of the sheath. Moving the slider forward in time, one observes that the electron profile predominantly oscillates in response to the wall potential, while the ion profile evolves considerably more slowly, particularly outside the sheath regions. The 1-D fields can be seen in Fig. 483.
image 1

Fig. 483 Plots of various 1-D field quantities showing the final state of the run.

Further Experiments

Now that we understand some of the basic physics of the discharge, we are in position to apply physics-based particle loading methods to hasten its eventual convergence to steady-state. The underlying principle here is to identify the ‘slow’ processes involved in the evolution of the discharge toward steady state, and then alter the loading to more closely mimic the state to which the plasma is being driven. While we cannot entirely predict the parameters of the steady-state, it is not difficult to at least get some idea of how the simulation is evolving and adjust the particle loads accordingly. We have already observed a number of physical processes of possible relevance:

  • initial electron loss and the establishment of ambipolarity
  • the slow decay of the total ion and electron population following the initial electron loss
  • the rapid response of electrons to applied electric fields, particularly in the sheath region
  • the slow evolution of ion density profiles.

Of these, we will primarily consider the ion profiles; the high mobility of the electrons suggests that electron profiles will adjust correspondingly on much shorter timescales. Additionally, since the strong electric fields in the plasma sheath region are screened out via Debye shielding as we move away from the walls, it seems clear that profile adjustments in the bulk plasma (where the driving electric fields are weakest) will ensue more slowly than in the plasma edge. We therefore concentrate our attention first on obtaining an approximately correct value for the ion density at the center of the domain.

From the ‘Data View’ menu at the top, select ‘1-D Fields’ again, and set the plot variable to ‘<None>’ in plots 2, 3, and 4. In Graph 1, set the plot variable to ‘He1Density’ and again move the timeslider on the bottom right of the window. Observe that the central ion density appears to be steadily but slowly rising; from its initial value of 5.12e14 1/m^3, it rises to 6.0e14 1/m^3 by the end of our comparatively short run. In addition, a rapid decrease in density near the walls (associated with the plasma sheath) has lowered the edge densities to around 1.2e14 1/m^3. The ion density can be seen in Fig. 484.

image 1

Fig. 484 The He ion density at the end of the run.

Iteration 2

Now we’ll create a new run which accounts for the central density buildup at the outset, by returning to the Welcome tab. The run we just did should appear in the list of Recent Simulations; we’ll want to click on that entry in the list and then hit the Copy Recent button on the bottom right. In the Select Files to be Copied window, the default settings are fine; just hit Select, choose a name and location for your new run, and save the files. We are automatically taken to the Setup Window, where we can adjust the density. Let’s try a value half again as big as the previous value; set PLASMA_DENS to 1.5 * 5.12e14 1/m^3 = 7.68e14 1/m^3 and Validate using the button at the upper right. Move to the Run Window and again adjust computation parameters (number of cores, etc.) as before before clicking Run.

When the simulation completes, return to the Analyze Window and rerun the ‘computePtclNumDensity.py’ script as before, inputting the new simulationName appropriate for the most recent run and again running the script for both the ‘electrons’ and ‘He1’ species. Then return to the Visualize Window. In the Data View menu, again choose 1-D Fields, set the plot variable to ‘<None>’ in plots 2, 3, and 4, and again select ‘He1Density’ in Graph 1. Moving the timeslider this time reveals that the central density continues to rise somewhat, but not as much as in the previous case. We could continue to raise PLASMA_DENS until this behavior ceases, but will forego that for the moment. We also observe that the edge density plummets sharply, ultimately nearing the 1.2e14 1/m^3 measurement we observed in the previous run. These observations suggest that a spatially dependent profile peaked at the center may be more appropriate. The ion density can be seen in Fig. 485.

image 1

Fig. 485 The He ion density at the end of the run.

Iteration 3

Let’s create yet another new run by returning to the Welcome tab and repeating the Copy Recent procedure we did a moment ago. When we are taken to the Setup Window, we now want to focus on the DENSITY_RAT parameter. Rather than loading the particles uniformly, we will now impose a symmetric, parabolic function (with value 1 at the domain center, and falling off to value DENSITY_RAT at the walls) as a probabilistic restriction of the particle load. Explicitly, this function has the form

\[f(x) = \frac{4x}{L}(1-\frac{x}{L}) + D(1-\frac{2x}{L})^2\]

where L = LENGTH (the extent of the spatial domain, as given in the input file) and D = DENSITY_RAT. The probability of generating a particle at point x is f(x); we thereby generate fewer particles near the walls. We’ll choose DENSITY_RAT as 0.1667 (i.e. 1/6) to begin with and see how things evolve. Make this change and save in the Setup Window. In the runtime options in the Run, again adjust the computational parameters (number of cores, etc.) as before, and run the new simulation.

When the simulation completes, we will need to repeat a number of steps we have done before:

  • Return to the Analyze Window and select the ‘computePtclNumDensity.py’ script, inputting the new simulationName appropriate for the most recent run and again running the script for both the ‘electrons’ and ‘He1’ species
  • Return to the Visualize Window, again choosing 1-D Fields in the Data View menu, setting the plot variable to ‘<None>’ in plots 2, 3, and 4, and selecting ‘He1Density’ in Graph 1.

As we move the timeslider, the evolution of the ion density profile is now much less pronounced, suggesting that our initial particle loading more closely resembles the steady-state profile. Yet if we switch the Data View menu at the top left to look at History, we can see that the particle densities have not yet come to steady-state. This is done by again repeating steps we have performed previously:

  • Set the plot variable to ‘<None>’ in plots 3 and 4
  • Set the plot variable to ‘physIons’ and ‘physElec’ in plots 1 and 2 respectively
  • Set the Location variable in plots 1 and 2 to the same value, e.g. ‘Window 1’.

Now we can see that both ion and electron densities decrease slightly as the simulation proceeds; a slight excess of plasma density is slowly being transported out of the system. This excess will persist until a balance is struck between the production rate and the rate of wall losses associated with a given profile; we are slightly over-predicting the density production and the plasma thus adjusts itself in such a way that more particles are expelled. The electron and ion populations can be seen in Fig. 486.

image 1

Fig. 486 The electron and ion populations versus time.

Lowering the PLASMA_DENS parameter to reduce the initial particle count is one way of dealing with this problem, but a more elegant solution which takes advantage of the more rapid physics processes occurring near the walls can also improve our initial particle loading. Recall that the initial behavior of the plasma, even before the ambipolar phase, centers on the establishment of the plasma sheath and the strong electric fields associated with the sheath. This behavior happens on a timescale much faster than any process in the ambipolar phase, and is characterized by high electron losses to the wall.

Iteration 4

Suppose that rather than loading particles all the way to the wall, we reduce the overall plasma density by leaving a gap near the wall where neither electrons nor ions are present. What happens in the discharge? The electrons, being highly mobile, rush to fill the gap, but rather than immediately being lost to the wall, they instead produce strong electric fields at the plasma edge which begin to modify the ion profile and bring about ambipolarity. If the gap is sufficiently large, the collisional production of ions and electrons will begin before appreciable wall losses ensue, and we can thus assess the relative rates of production and loss fairly early in the simulation. Since the electrons are highly mobile, let’s treat the average electron population as a measure of how well we’ve achieved this balance; net electron production as we move into the ambipolar phase means that our gap is too low (we have removed too much density), while net losses mean that our gap is insufficiently wide. As the profile shapes near the walls tend to adjust themselves fairly quickly (due to the larger electric fields in this region), we can in this manner obtain approximately correct values for the total ion and electron populations at the simulation outset.

Return to the Welcome tab and create another new run, using the run we just carried out as a basis. In the Setup Window, the parameter ‘SHEATH_GAP’ describes the physical width of the particle-deficient region we’ll introduce near the walls. We’ll need to do a bit of work since our profile is now parabolic. Recall that the probability function we’re using for the parabolic load has the form

\[f(x) = \frac{4x}{L}(1-\frac{x}{L}) + D(1-\frac{2x}{L})^2\]

where L=LENGTH and D=DENSITY_RAT. If we average f(x) across the simulation domain, we obtain

\[<f(x)> = \int_{0}^{L} \frac{f(x)}{L} dx = 1 - \frac{(1-D)}{3}\]

The function <f(x)> is proportional to the average species population; introducing a gap of width S = SHEATH_GAP to reduce this average population modifies this calculation such that

\[<f(x)> = \int_{S}^{L-S} \frac{f(x)}{L} dx = (1 - \frac{2S}{L}) - (1 - \frac{2S}{L})^3 \frac{(1-D)}{3}\]

Assuming we want the new <f(x)> to equal some fraction B of the original population, we can find the appropriate sheath gap by numerically solving the cubic equation

\[B(1 - \frac{(1-D)}{3}) = (1 - \frac{2S}{L}) - (1 - \frac{2S}{L})^3 \frac{(1-D)}{3}\]

for the SHEATH_GAP parameter S. (Alternatively, we could just make a guess for a sensible value and adjust accordingly.) For our case, with D = 1/6; letting B = 3/4 yields three real solutions to the cubic, but one is negative (and thus unphysical) while another is greater than L (and thus also unphysical). The remaining solution yields a value 0.0133 for SHEATH_GAP, so let’s try that as a starting value.

  • Run the simulation with the modified SHEATH_GAP parameter, with other parameters the same as were used previously
  • Return to the Analyze Window and select the ‘computePtclNumDensity.py’ script, inputting the new simulationName appropriate for the most recent run and again running the script for both the ‘electrons’ and ‘He1’ species
  • Return to the Visualize Window, again choosing History in the Data View menu, setting the plot variable to ‘<None>’ in plots 3, and 4, and selecting ‘physIons’ and ‘physElec’ in plots 1 and 2 while giving the Location variable the same value.

The electron and ion populations can be seen in Fig. 487.

image 1

Fig. 487 The electron and ion populations versus time.

The average electron population, after the initial wall losses, remains relatively constant, so we have struck a good balance between ionization (an electron source) and wall loss (an electron sink). Modified loading techniques have enabled us to observe and balance the rapid physics processes which ensue as electrons rush to fill our imposed sheath gap S. Such balancing hastens the ultimate convergence of the simulation to steady-state.