Utilizing VSim for Advanced Simulations in Plasma Physics and Design Engineering With Joe Theis

[00:00:00] Hello and welcome back to TWSS 2023 on our second day. The next talk that we have is a graduate student who is studying under Dr. John Cary right now at the University of Colorado. This is Joe Theis.
Great. Thanks, Colleen. Yeah. So I’m Joe Theis I’m currently at CU Boulder and today I’m going to give a presentation about PIC simulations of magnetron discharge.

Yeah, I’m going to start by outlining some of the applications of magnetron discharge, specifically sputtering. I’m going to go over some of the essential physics to this talk. I’ll introduce my research question, go over my simulations. They should set up and be some how I’m simulating this problem and then discuss some shortcuts to make this problem more computationally tractable.

And then finally, I’ll give you some results. and end with some conclusions and unknowns. Magnetron discharge is primarily used for sputtering [00:01:00] and so that’s to sputter coat thin films. In the image on the left, you can see a planar magnetron discharge where you see the characteristic glow of the plasma there in the ring and so that’s used to deposit thin films on various materials.

It’s used in metallization and semiconductor fabrication as well as optical coatings for things like photovoltaics shown on the right. And also even architectural glass, so they can span really large areas to really small areas, and it’s the most common technique used to deposit thin films. Here are some of the essential physics for this talk.

On the right, I’m showing a cross section of a planar magnetron discharge chamber, and so it’s a cylindrical chamber, and so this is the symmetry axis in the middle. And I’m going to start by introducing DC diode sputtering, and so that’s without the magnetic field, and that should motivate the need to introduce these magnets.

[00:02:00] So here in this image, we have a cathode on the bottom and an anode on the top, where we negatively bias the cathode, which gives us an electric field. And when this electric field is sufficient, we can achieve a glow discharge. And that’s what we’re looking for. Or we have this self sustaining cycle where electrons get accelerated and ionized, our background gas, so this chamber is filled with some neutral gas, in our case it’ll be argon.

And then those argon ions will be accelerated by the electric field where they strike the cathode and can produce more secondary electrons. And that’s how it self sustains. But what makes this process useful for sputtering is that those ions also emit sputtered atoms from that cathode. So that cathode is made of some target material like titanium or copper.

And so the ion comes in, strikes that material, frees up that copper atom, and then it travels holistically where it lands on this substrate. And we start growing a thin film up here. And so that’s the goal of it. And so that can be accomplished in diode sputtering. However. To do that, you [00:03:00] need high pressures, and that’s because the electrons need to ionize a sufficient number of atoms before they cross our domain.

And so due to those high pressures, these sputtered atoms have a high probability of collision, so that decreases your deposition rate. So you can’t grow these films as fast, and they’re also less uniform, and you get a lot of sputtering on the other parts of the chamber where you don’t want it. And so due to those limiting factors of this DC diode sputtering, non sputtering was developed, and so that just involved adding these permanent magnets to the back of the cathode, and you can see that produces this B field as I’ve shown in the figure here.

And so what this allows us to do is trap the electrons. So electrons will be doing an E cross B drift as muthily in this image here. And so that prolongs their lifetime, meaning that they can execute many mean free paths at much lower pressures. So now we can lower the pressure by an order of magnitude.

And again, that gives us a much higher deposition rate because now those sputtered atoms can just travel without any collisions and deposit on our substrate. [00:04:00] Yeah. And so we get a higher deposition rate by an order of magnitude. Dude, and it also operates at lower pressure or lower voltages, so it’s more efficient.

That process works well, but as with everything, we’re always trying to improve it and motivated by the image on the right, that’s a photonic integrated circuit that I’m holding. And you can see the level of detail in that image on this integrated circuit. So the demands for all of these semiconductor fabrication techniques are getting tighter and tighter as we want more smaller details on these devices.

And so therefore we want to improve these different devices involved in chip manufacturing. I’ve broken out some of the technological goals. Two categories. One is efficiency manufacturing. The other is film quality. In efficiency of manufacturing, what is it? It’s improving the deposition rate, decreasing the wall power, increasing our target utilization.

How much of that cathode block of material do we actually use? And increasing the sputter yield. This is important, but reduces energy use and cost of [00:05:00] fabrication. You can think of that as being especially important in optical coatings that cover really large areas. How can these be accomplished? You could optimize the B field to increase electron trapping with the goal of decreasing the voltage.

Another broad The quality of the films, the film quality includes things like its thickness, uniformity, density, and purity. This is important, again, for the different integrated circuits that are being printed. And how people are attempting to address this now is using ionized sputtered atoms. And so this is extending the regime of DC magnetron sputtering, where now you’re applying a pulse.

This pulse actually does a better job of ionizing those sputters, which creates a denser film. And so that’s an eventual goal of my research is to model those which are really high density and high currents. Yeah. But for now, I’m interested in investigating the IV characteristic of these magnetron discharges.

As I mentioned before, the voltage is critical for our controlling these sputters ’cause that controls your [00:06:00] current and therefore the sputter rate and your growth of your film in these chambers. And right now we only have experimental results and people. I’ll fit different curves to these models, but there’s no complete quantitative model explaining this IV curve for magnetron discharge.

The IV curve also expands three orders of magnitude in current, so there’s lots of different regimes or different physics are going to be important. And in terms of simulation, there’s no converged PIC simulations that cover any, yeah, any range of this IV curve. And so that leads me to my research question, which is, can we quickly simulate industrially relevant magnetron discharges with PIC?

And how I propose to answer this is by benchmarking our PIC simulation against an experimental IV curve. It will address the question on the top, because if we can accurately benchmark, give that question an affirmative yes, it’ll also help us identify some relevant physics in the problem. We want to explore other regimes as well.[00:07:00]

And so that is the project I’ve currently been working on. And so here’s the simulation set up to simulate these magnetron discharges. I’m simulating in 2DRZ cylindrical. And so now the Z symmetry axis is horizontal, and then we have the radial direction going vertically. Everything inside this box is the domain that I’m simulating.

And so you can see that I have a magnetic field, and so that’s generated by a magnetostatic field solve. Before my simulation, I solve for this field, and then that’s a constant throughout. Then, to introduce the electric field, I bias this cathode, which is on the left. And that’s biased based on this external circuit where I have this external resistor in series with this voltage short so that the voltage applied to this cathode, it’s just that voltage minus the current collected at the cathode times this resistance, and that gives me a feedback mechanism to control the voltage, which I’ll talk more about later.

I filled this with a uniform neutral argon gas, and then with that argon gas. [00:08:00] So I’m considering electron gas collisions, elastic ionization and excitation at the cathode. I have ion induced secondary electron emission and electron reflection and the other boundaries are absorbing and I’ll talk a bit more about those boundary conditions.

This anode is the ground and then this cathode is negatively biased. So a bit more about the external circuit in this image on the right, I have an I. V. Characteristic with current on the Y axis and voltage on the X axis. And then in red, I’m trying to show some unknown I. V. Curve for the magnetron that I’m trying to discover.

And so that could have regions where there’s negative differential resistance, like I’m showing here in this regime, and regions where there’s positive differential resistance up here. And so that’s what I’m trying to develop with my simulations, and how I do that, how I access these different points is with this external circuit.

And that’s what I’m showing in blue. So these blue lines represent the IV curve for just a resistor in parallel with that [00:09:00] voltage source. So that’s going to be a straight line where the slope is equal to 1 over that external resistance. And so that forces the simulation to land on this regime in the IV curve.

This region is also unstable, so this external circuit stabilizes this regime. And then by decreasing the resistance between different simulations, I can actually sweep along this curve so that I can populate it and generate this IV curve.

So now a bit about the simulation scales. For a 0. 3 amp discharge the plasma density is about 10 to the 18 per meter cubed. The temperature for the electrons is about 6 eV. The magnetic field is 400 Gauss. The background pressure is 5, 5 millitor. And a typical domain size is 3 by 3 centimeters for this RZ cross section.

And we reach steady state in about 10 microseconds. So I’m giving you some length scales in this top table normalized by the Debye length, which is the smallest length scale. And [00:10:00] that’s what sets my cell size. So you can see that the electron gyro radius is 10 times that Debye length at this density.

And then the ion gyro radius is much bigger. In fact, it’s bigger than the device. And again, that’s why the ions are considered unmagnetized. And then for time scales on these devices, the smallest time scale is the electron crossing time of this Debye length. And so I’ve normalized all the other timescales based on that, since this is my time step for these simulations.

You can see that this density, the plasma period is 62 times that crossing time. Larmor period is 500 times, and the ion crossing is 400 times. And this kind of represents a response time of the plasma to different inputs at the cathode. I’ve also given a lot of the collision periods for this. Vice on the right side of this table, and this is just a reminder that the mean free paths are really long for these collisions, and so they’re really rare, but they’re important because they drive all the processes that are going on inside this magnetron discharge.[00:11:00]

You can see that we have all these electron collisions, which we are modeling at 2000 to 17, 000 electron crossings. And then I’ve also included the ion collision periods here, and you can see those are an order of magnitude higher, which is what I’ll later use to justify ignoring those. The last column on this table is the steady state time.

So this is the number of time steps it takes to get to this 10 microsecond mark in physical time, and that’s 5 million time steps. So with the performance I was seeing with just the default decomposition and this simulation setup being 3×3, It would take about 20 days on 2000 cores. And so I feel like a bit of a whiner after Faircloth’s talk in which he was talking about 16 week or 13 week simulations, but this still, if I want to sweep along the entire IV curve, that would take a lot of time.

That brings me to some shortcuts to an answer, because some of these shortcuts certainly introduce error, and I’ll talk about that. Some of the low hanging fruit [00:12:00] is I ignore sputtering and the deposition of those sputtered atoms on the anode. In the regime that I’m looking at, which is relatively low currents and materials with low sputter yields like titanium, this is, this simulate, this approximation has a negligible effect.

And that’s because rarefaction, which is a process by which those sputtered atoms deplete your neutral background gas, that’s going to be relatively low in this regime. Also, the ionization fraction of those sputtered atoms is very low for this regime as well, since we’re not really in that high PIMS, highly ionized regime that I talked about earlier.

Second approximation I’m making is I’m ignoring those ion collisions, like I mentioned on the last slide, they’re very rare. This does have an effect. It decreases the plasma density by about 10%. And in this table on the bottom for the Monte Carlo framework, I show that dropping, going from 3 collisions to 2 collisions to 1 collision, you get about a 10 percent speedup for each one of those drops.

Those speedups are more significant for the [00:13:00] reactions framework. And so I drop the ion collisions. Then in terms of the domain, the main goal is to reduce the number of cells. And therefore the number of particles as well. The first thing I do is I increase dr by a factor of 2 compared to dz. This is commonly done in PIC simulations of magnetron discharge.

It’s motivated by a couple of factors, but clearly strong variation in the fields occurs along Z in this sheet. It does increase the plasma density again by about 10%, but it results in about a three times speed up in the simulation. And again, and that’s due to many fewer particles due to this smaller domain and the fact that we can increase DT because DT is limited by DR for these 2D RZ simulations.

Then the other two things have to do with actually shrinking the domain, which I’m showing on the right. On the right is the full three by three centimeter domain where I show the potential and then I show the plasma density here. And so what I do is I take this full three by three centimeter and I bring it to a 20 millimeter by [00:14:00] 15 millimeter domain.

And I put a Neumann boundary condition on this upper axis for the potential. And you can see that looks reasonable if you look at the potential in this unperturbed field. And also there’s very few particles at that high boundary. Then at this upper boundary, I bring it in because the potential is constant at that point.

And again, there are a few particles. This is logarithmic scaling for the density plot over here. And also this 15 millimeter gap distance is within the range of experimental. Devices as well. So didn’t feel too extreme to do that. Yeah. The next strategy for speeding these things up is looking at the step response of this system to the turning on of the voltage to basically ignite this discharge.

And so when looking at the step response of any kind of device. It’s all about balancing the rise time and the overshoot and ringing. So this plot in blue on the right, this shows a constant current simulation where my power source isn’t operating in a constant current mode. And this shows a lot of the good features.

So the [00:15:00] rise time basically represents the time between when you turn it on. When it reaches that steady state behavior. And so that’s this region. And you can see that this blue curve has the fastest rise time on this plot. But with that fast rise time comes a lot of overshoot, which is the peak of this extending beyond the steady state number of ions that we expect.

And then this overshoot leads to ringing. And so that ringing results in this oscillation around your steady state. And so the ideal step response is really fast rise time with no overshoot and ringing, but it’s always a balance. And so the three different techniques I’ve shown in this plot are the constant current, which I already outlined.

Then there’s a constant voltage source. That’s shown in green here, and these are different simulations, so you don’t have to worry that they’re giving different results because they don’t. But this constant voltage source has a really long rise time, but then no overshoot and ringing. And then this finite resistance, which is the external resistance model that I’ve already discussed.

It’s a happy me. [00:16:00] Yeah, more. The rise time is intermediate, and then it has a little overshoot ringing, but it’s much smaller than the constant current source. And then finally, one thing that I’ve investigated is using a time dependent resistance, and this can allow you to get that fast rise time and low ringing.

And so that’s a technique of used as well, and I’ll show some results. So that later on. And one more thing I want to mention about the external circuit is that constant voltage, you can get to an answer, but in some cases where you have a negative differential resistance in the IV curve, it destabilizes the discharge since there’s no feedback.

So that’s what I’m showing in this block, this red line to the left of it. This is time to the left of it. There was an external resistance. And so that was providing the stabilizing feedback. And there was this nice steady state reach for the number of ions. But then to the right, I restarted the simulation with a constant voltage at whatever the steady state voltage was.

And you can see that the number of ions starts steadily increasing, and the simulation goes unstable. And in some cases on the IV curve, you need the external [00:17:00] circuit in order to stabilize the discharge. All right, this is the last technique I had to, and this, in VSim, you’re able to have a lot of freedom when you decompose your simulation into different grids.

And so that allowed for this development of the recursive decomposition algorithm. And what I’m showing on the left is a slice decomposition where the slices have to go all the way through the domain. And what’s different about the recursive decomp is that those slices can end and you can have more freedom to.

Doing these domains to balance the load in the different MPI ranks for whatever simulation you’re running. I know Daniel Murray talked about this yesterday, but I thought I’d touch on it as well too. So this is showing the particle density in the background. And you can see how the recursive decomposition does a better job of mapping to that particle density.

The bottom stat on this plot shows the maximum and minimum relative work in a domain. For all of these domains, for each of these different decompositions. [00:18:00] And they’re relative to the average work in a domain. So in the ideal decomposition, all these values would be one and you can. See that the values in the recursive decomp are much closer to one than in the slice decomposition, and I see those results bear out in terms of performance.

So here is a table with a lot of data in it, but I think what we only need to select out here is the full speed up that I got. Originally, I was using a default decomposition, and the simulation would normally be run on 32 cores. This table is showing the strong scaling. But for the performance, which throughout this talk will be given in core microseconds per particle step, that value is 6.

4 for the reactions and the default decomposition. And by moving Monte Carlo and using a recursive decomposition, I was able to drop that to 1. 8. So that’s a three times speed up. with this recursive decomp. And this also isn’t an ideal decomposition, since I was assuming that the solve, or that the simulation was particularly, that [00:19:00] the work was dominated by particles.

I was using this value of a equals four, and that value comes into this, basically how we load balance. We have some approximation where we say the work per cell is given by one plus a times the number of particles per cell. So this is basically a ratio of the work done between the field solve and the particle push.

And so what I did is I varied the number of particles per cell, and then I looked at the core microseconds per cell step. And sure enough, it was nice and linear, and therefore I could fit it to this line. And from there, I get a equals point eight, which is different than what Daniel Murray presented yesterday, which was the value of equals four.

So in my simulation, this is definitely gonna be simulation dependent. It’s dominated by the field cell. And so with this value of a equals 0. 8, I get even a more significant speed up for the recursive decomposition.

So finally, I’ll go over some [00:20:00] results for a 0. 3 amp discharge. I guess I want to start out by emphasizing that earlier I quoted computational time of 20 days on 2000 cores. And so I was able to run this simulation to steady state in about six days on about 500 cores. So that’s about a 16 times speed up, which is very significant and makes this simulation a lot more tractable.

But this is a 0. 3 amp discharge, which is right in the middle of the regime I want to be looking at. This is definitely an industrially relevant magnetron discharge. And this plot on the bottom left shows Voltage of the cathode as a function of time. I restarted this from a lower current. That’s a technique I’ve been using as well.

So I restart from a steady state achieved at a lower current for my next simulation. And then I do a time dependent resistance, which basically ramps that resistance down to the next resistance I want to simulate. If you remember my external resistance plot from before. And then I allow it to evolve freely at that constant resistance, in [00:21:00] which case I get a little bit of overshoot, but the ringing settles down really quickly before I get to steady state voltage.

and it doesn’t look totally clear that it’s a steady state on this voltage plot. But if you look at the number of ions, it’s a bit more clear that there’s this rise time and then some ringing. But we’re in the steady state regime. And so this is just one of the many simulations I’ve been running, but I thought I’d just show some basics of what the different fields look like as well.

On the left, I’m showing the magnetic field. And so we have a magnet in the middle, that’s a cylinder, and then we have this annular magnet on the outside, and you can see that by the peaks and the magnetic field intensity given by the color plot. Magnetrons are characterized, or the magnetic field in these magnetrons is characterized by the maximum tangential field along the cathode.

And in this case, that’s 350 gauss for the simulation. The potential is shown in the middle where you can see Do that really thin sheath. And the cathode is at 280 volts here. And then finally we have the plasma density. This is linear scaling, [00:22:00] but we get this real peak. And this is called the ionization region in the middle here, where the majority of ionization occurs.

And then the ions fall from there into our cathode.

And so finally, I just showed a 0. 3 amp simulation. And so that corresponds to this orange point on the far right of this curve. I guess after getting you guys used to having IV curves with current on the y axis, I’ve decided to switch it to keep it interesting, but all of these points represent independent simulations.

Sometimes I restart them to achieve the higher point, and so this is my data in from PIC in orange, and then I’m showing some experimental results with a titanium target in blue. They look very different. But I think there’s a lot of sensitivity to magnetic field and different effects. And also there’s really low resolution in this experimental curve on the left.

Clearly I have this big feature here where I have this peak in voltage, and then it settles down and I get this [00:23:00] flat section. All of these results are on the peaks. Relatively well converged. And so I’m currently working on explaining the physics behind why this peak in the voltage occurs. I’ve included this plot on the right, which is of the IV curve of a DC glow discharge, no magnet.

So this isn’t the yeah, so this isn’t a magnetron discharge ID. But I think it looks similar to what I’m seeing in so far with my pig simulations with this first hysteresis where it rises up. And then it settles down to this flat section. So I’m curious if I’m discovering a similar hysteresis process for magnetron discharge, which hasn’t been reported.

And so yet to be seen is whether I get this eventual rise in the voltage, which I expect to see at higher voltages. And yeah, I’m hoping to see that in the future.

Yeah, some conclusions from this work is that we can simulate industrially relevant magnetrons in a reasonable [00:24:00] amount of time. Power sources must be informed by the I. V. You have to know where on the I. V. You’re intersecting based on your power source and whether or not you need to stabilize it. The shortcuts can lead from 10 to 20 times speed ups that I’ve mentioned before, but some of them come different errors that need to be quantified.

The best shortcuts include the recursive decomposition shrinking the domain size. And managing the plasma response would be external circuit. Some unknowns. So first of all, the low current results are converged. As I mentioned, high current results don’t appear to be converged in terms of the number of particles per cell.

And I think it remains an open question of if these points will converge. This plot on the right shows a recent result from boss at all in 2022, where they were simulating a one dimensional RF discharge. And what they saw this end is the. Number of particles in the domain. So more particles moves us to the left.

So as they increase the number of particles in the domain, they actually saw a divergence of the result [00:25:00] here, given as the plasma density and so far, I’m seeing a similar concept, whereas I increase the number of particles. I get this divergence, which is very concerning. But beyond that, I want to look at physics explanations for the characteristic and get a solid explanation for that.

Yeah, I think I have some more results, but I think I’ll leave it there and open it up to questions. Thank you. And this work was done with Gregory Werner, Thomas Jenkins, Daniel Mahan, and John R. Carey.

Thank you so very much, Joe. That was great. If anyone has questions, we have time for one or two. You can type them in the Q& A box.

At the bottom of your screen, and we can read them aloud.

And I will say a few things while we’re letting people type in their questions. We will be making videos of all of these presentations. And they will make their way to the TechX YouTube channel. It takes a little bit of time. It might be one to two weeks before they are all up and available. But we will have them there for [00:26:00] you to go back.

Another thing I want to mention is if you have not evaluated VSim, And 1 of our latest versions, we are up to be some 12. 2. And so if you are a current customer with an older version, or if you haven’t evaluated us in the last year go online, request an evaluation and let’s get you playing around with the latest and greatest version that we have for you

1. I should say that latest version. It looks like we’ve got 1 question. And it starts out with a compliment to you. Great talk, Joe. Do you observe grid heating in the radial direction when you reduce the radial grid resolution?

Thanks, Luke. And yeah, I do see some grid heating. And again, that results, I talked about it in the talk briefly, but I see about a 10 percent increase in the plasma density.

And I think that increase in plasma density is due to some form of grid heating. Due to that radial [00:27:00] direction, which leads to more ionization, therefore a higher plasma density, and yeah, now that I’m struggling with convergence at these high currents, that’s something I’d want to check again and see if I can really play the games with the radial direction like that.

But to clarify… There’s, you’re not seeing an exponential instability in terms of the plasma density? Oh, in terms of the voltage, it does appear to be exponent, exponential right now. Like it, remember how it’s diverging, but I haven’t pushed too high. This is just, yeah. So I could certainly get more particles per cell.

But I think, I’d definitely like to hear from different Tech X people who have modeled CCPs, whether or not they see convergence based on the particle per cell count.

All right, I think

that’s all the questions that we’ve We have, we are going to take a short break. We will be back at 11 a. m. Mountain time for our next talk [00:28:00] by Dr. Daniel.