NetLogo User Community Models

(back to the NetLogo User Community Models)

Download If clicking does not initiate a download, try right clicking or control clicking and choosing "Save" or "Download".(The run link is disabled for this model because it was made in a version prior to NetLogo 6.0, which NetLogo Web requires.)

ABSTRACT

NetLogo has been used here to solve a typical numerical problem, namely a non-linear differential equation. In this case the interacting agents are simply the solutions of the equation at previous time steps. The equation considered here describes the behavior of a superconducting Josephson junction, however it should be simple to adapt the program to different differential equations.
A full report of the present project is in preparation and paper and web-based references to it will be added to the documentation as soon as possible [Maggi03].

WHAT IS IT?

A Josephson junction is a superconducting electronic device based on the tunneling of Cooper pairs across two weakly coupled superconductors. Besides its current practical applications, the Josephson junction is important from a physical point of view because it has been the first device showing quantum mechanical effects on a macroscopic scale.

According to quantum mechanics, a single electron has a small but finite probability to tunnel through a potential barrier higher than its kinetic energy. In a Josephson junction, however, also pairs of electrons ("Cooper" pairs) can tunnel through the barrier, usually made by a thin oxide barrier (1-2 nm) sandwiched between two superconducting films. The physical reason is that the two superconductors are so close together that the wavefunctions of the Cooper pairs of the films overlap and the system behaves as a single "weak" superconductor.

What is observed experimentally is the flow of a current across the junction with no voltage drop (V = 0), at least until the current reaches a "critical" value Ic. When the current is further increased, also single electrons originated by the breaking up of the Cooper pairs begin to tunnel the barrier. Therefore, the potential difference V between the two superconducting films becomes ­ 0 and a state is reached where the junction behaves as a resistance, Rn.
Further details about the basics of the Josephson effect at an undergraduate level can be found in [Langenberg66], [Anderson70].

The behavior of a Josephson junction can be modeled by a simple equivalent circuit consisting of the parallel connection of a resistance R, a capacitance C and a Josephson supercurrent Ic sin(phi), where phi is the phase difference of the wavefunctions of the Cooper pairs across the junction. The circuit is assumed to be current biased (using a dc and an rf bias) and the voltage V across the junction is proportional to the time derivative of the phase phi. The circuit resistance R is the normal state resistance Rn while C is the capacitance of the two superconducting electrodes separated by the thin oxide barrier.

This model, known as the resistively shunted Josephson (RSJ) junction model, is accurate enough to gain a basic understanding of the junction behavior. Despite its simplicity, the empirical RSJ model is the limiting case of the complete quantum mechanical treatment of the Josephson effect. The most important aspect left out of the model is the spatial variation of the phase, which is related to the effect of a magnetic field on the junction.
The RSJ model is also able to show more advanced effects such as the chaotic phenomena observed in Josephson junctions.

The differential equation representing the RSJ model, using the McCumber dimensionless parameters [McCumber68], is

d^2 phi d phi
beta-c ------- + ----- + sin(phi) = alpha-dc + alpha-rf sin(omega-rf tau) (1)
d tau^2 d tau

where phi is the phase difference across the junction, tau is the normalized time, beta-c is the damping parameter, alpha-dc and alpha-rf are the dc and rf current bias in units of Ic, and omega-rf is the normalized rf frequency. Eq. (1) is also the equation of motion of a damped pendulum driven by a constant torque plus an oscillatory sinusoidal torque. In this case phi represents the angle of oscillation from the vertical.

A higher-order differential equation can always be transformed in a set of first-order coupled differential equations. In particular, (1) can be reduced to

d y1
----- = y2 (2a)
d tau

d y2 1
----- = ------ [alpha-dc + alpha-rf sin (omega-rf tau) - y2 - sin(y1)] (2b)
d tau beta-c

where y1 = phi and y2 = (d phi / d tau). These equations can easily be solved using standard numerical methods.

Three simple algorithms have been used to solve Eqs. (2). The first is the basic Euler's method, well known for its simplicity but also for its lack of precision. The second is the second-order Runge-Kutta method: this method has been preferred to the more powerful fourth-order Runge-Kutta method because it is much simpler to implement in NetLogo and gives results precise enough for the present purposes. Both methods have been implemented in the "individual" model described in the next section. Finally, a finite difference method has been used in the "cooperative" model also described below.

HOW IT WORKS

To show that it is possible to use different approaches to write a NetLogo program, two alternative NetLogo models have been developed. In both cases, the Graphics window shows the time dependence of the derivative of the phase, phidot, which is the physically interesting quantity, being proportional to the voltage across the junction. In addition, a phase space plot, i.e. the dependence of phidot on the phase phi, can be shown (at the expense of simulation speed).

The "individual" model shows the approach one would take using a standard programming language such as C, Fortran or Basic. The simulation starts with a single turtle. For each simulation step, the equations (2) are solved using either Euler or Runge-Kutta methods, and a new turtle is created whose (x, y) coordinates are proportional to the current simulation time, tau, and to the current value of phidot, respectively. To avoid to have too many living turtles, which would slow down the simulation and clutter the Graphics window, the turtles age and, after a finite time of living, die.

The "cooperative" model is closer to the "proper" NetLogo approach to programming. In this case the simulation starts with a population of turtles equal to the width of the Graphics window and evenly distributed along the y = 0 axis. The x position of each turtle is kept fixed, while its y position is proportional to phidot, i.e., to the voltage across the junction, calculated using a finite-difference method. With this approach, no new turtles are created during the course of the simulation nor aging turtles die. This leads to a simpler code and speeds considerably up the simulation with respect to the "individual" model.
In reality, this approach is weakly "cooperative", because the single turtles interact only because the vertical position of the (n+1)-th turtle is calculated knowing that of the n-th turtle. An alternative, truly cooperative approach for this model would be very interesting, and contributions are welcome.

HOW TO USE IT

The NetLogo model interface is divided into four main areas: a set of sliders on the left change the simulation setting, while several buttons, sliders and switches on the bottom control how the simulation is run. The Graphics window shows the main result of the simulation together with the plot of the phase-space on the right. In the following all the variables will be in uppercase, to make them stand out of the main text; however in the interface window the variables are in normal lowercase.

The top three sliders set the simulation settings: the simulation time step DTAU and the values of the initial conditions for the solution of the differential equation, PHI-0 and PHIDOT-0. The four sliders below set the main junction parameters of eq. (1), dc current bias ALPHA-DC, damping parameter BETA-C, rf current bias ALPHA-RF and normalized rf frequency OMEGA-RF.

The button SETUP sets up the NetLogo model by creating a population of turtles equal to the width of the Graphics window and distributing them on each patch with y = 0. The procedure also clears and sets up the phase-space plot (see below). The button DEFAULT-PARAMS sets the default simulation parameters, useful to start experimenting with the model.

In the "individual" model, the switch RUNGE-KUTTA? selects the algorithm used to solve the differential equation: when ON, the second-order Runge-Kutta method is used, when OFF the simpler Euler's method is selected. This switch is not present in the "cooperative" model, which uses only a finite-difference method to solve the differential equation.

The button START/STOP is a forever button used to start and stop the simulation, while the button SINGLE STEP advances the simulation by one time step. The SLOWDOWN slider is used, when necessary, to slow down the simulation in order to better observe the time dipendence of the waveform.

The Graphics window shows the dependence of PHIDOT, namely the voltage across the junction, on time TAU. Each dot in this window represents a time step and hence leads to the artifact that the waveform widens or shrinks by changing the time step DTAU. In this respect, the Graphics window works more or less like an oscilloscope (sorry, no trigger yet!) To improve the visualization of the waveform, a slider Y-AMPL has been added, that amplifies by the corresponding factor the vertical scale of the waveform.

On the right of the model window, the phase-space plot of the simulation is shown. This is a plot of the state variables phidot vs. phi and gives a very useful insight of how they are related during the course of the simulation. Since the plot slows down the simulation considerably, it is off by default and can be activated by means of the PLOT? switch. Also, a RESET PLOT button clears the plot and resets the PLOT? switch to its off state.

To use the model, one must first set the simulation settings and the junction parameters and then click on the SETUP and START/STOP buttons. The variation of one of the parameters during the simulation has an immediate effect on the calculated waveforms except, obviously, for the values of the initial conditions which affect the model only at the beginning of the simulation.

THINGS TO NOTICE AND TRY

Start the model with the default settings: PHI-0 and PHIDOT-0 = 0 and the time interval DTAU = 0.1. The parameters ALPHA-DC, ALPHA-RF and OMEGA-RF are all zero while BETA-C is set to its minimum value of 0.1. The waveform vertical amplification is 5x. In these conditions the junction has very little damping and the external current is zero. The value of SLOWDOWN depends on the speed of the computer used to run the simulation. It is better to set it initially to zero and increase its value until the waveform changes slowly enough to be clearly seen on the screen.

Press SETUP: a yellow horizontal line should appear in the Graphics window, showing the turtles evenly distributed along the y = 0 axis. Press START/STOP: nothing seems to appear. The normalized voltage across the junction is zero and therefore all the turtles remain at their position on the y = 0 axis.

Try to increase slowly ALPHA-DC, up to 1.0. No change in the turtle position is visible. A current of Cooper pairs is now tunneling across the junction, but the voltage remains zero because the junction is acting like a single superconductor and the current can flow with no electrical resistance and hence no voltage drop.

Increase further ALPHA-DC just above 1.0. Small voltage spikes appear, moving relatively fast towards one side of the Graphics window, in a manner similar to an oscilloscope with a slightly incorrect trigger level. Single electrons have started to tunnel across the junction barrier along with the Cooper pairs. For them, the junction is a nonlinear resistance, and hence a finite time-dependent voltage appears across the junction electrodes.
Changing DTAU, the waveform can be made almost still, like adjusting the trigger level on the oscilloscope. The exact value of DTAU depends on the speed of the computer used to run the simulation.

Increase ALPHA-DC further. The frequency of the spikes increases and their form becomes closer and closer to that of a sine wave. Likewise, also the level of the waveform increases, therefore the average voltage across the junction is increasing. Play with DTAU (and Y-AMPL) to experiment how the image of the waveform changes by changing these parameters.

Reduce ALPHA-DC to 1.2. Set OMEGA-RF = 0.2 and increase ALPHA-RF slowly from 0.0 to 2.0 observing how the waveform changes because of the applied rf current. It is better to start with a relatively large value of DTAU, around 0.1. On fast machines, set SLOWDOWN to a value that allows to see clearly the resulting waveform (for example, on an 800 MHz Macintosh G3, SLOWDOWN = 0.4 - 0.6 with the phase-space plot off is adequate).
The voltage spikes, evenly spaced at ALPHA-RF = 0, seem to group and the number of spikes in each group is related to the amplitude of the rf bias, ALPHA-RF. Play with ALPHA-DC, ALPHA-RF and OMEGA-RF to observe what happens by changing these parameters.

Finally, try **exactly** these parameters: ALPHA-DC = 2.1, BETA-C = 0.1, ALPHA-RF = 2.0, and OMEGA-RF = 1.5.
Also, set DTAU = 0.1 and Y-AMPL between 5x and 10x. The waveform is periodic, as is apparent from the phase-space plot. However, it is sufficient to increase ALPHA-DC to 2.2 and above to observe a completely different situation: the resulting waveform loses its periodicity and starts exhibiting chaos.
This is more apparent by looking at the phase-plane plot. In all the above examples, the phase-plane plots at subsequent periods of the rf bias overlap as a result of the existence of a periodic solution, while in this latter case the plot continues to change because of the chaoticity of the solution, ending up to completely filling a large region of the phase-space.
It is not possible to delve more deeply into the subject of chaos, here. The reader interested to a simple introduction to chaos in Josephson junctions can refer to [Kautz93]

EXTENDING THE MODEL

Among the many things that could be made to extend this model, it would be simple to add a plot of the Poincare' section of the phase-space, i.e., the values of phi and phidot at the beginning of each cycle of the rf drive. This plot provides unambiguous information about the form of the solution, whether periodic or chaotic [Kautz93].

Another trivial modification to the model would be to change the phase-space plot, using cos(phi) vs. phidot instead of phi vs. phidot, maybe by implementing a switch that sets up properly the scales of the plot, the axis labels, etc. This provides an alternative view of the phase-space plane, with the added important benefit that a periodic solution appears in this plane as a closed overlapping curve while a chaotic solution is easily recognizable because the resulting plot fills completely a region of the phase-plane.

NETLOGO FEATURES

NetLogo seems very useful not only for a problem with interacting agents but also for a purely numerical problem like this. The inherently parallel approach of the NetLogo language leads to a clean and simple to understand code. Of course, for heavy numerical work, a conventional language such as Fortran or C is still preferrable, for its raw speed and also because it is possible to run the simulation in unattended tasks, possibly collecting and analyzing automatically the results.

The real benefits of NetLogo are different: the language is very simple but also powerful and a complete simulation, comprising a fully interactive user interface, can be built with very few lines of code, making it an invaluable didactic tool or even a tool for rapid prototyping of an algorithm. The present simulation has been written with about 70 lines of code (not considering comments, blank lines and isolated square parentheses.)

Last but not least, NetLogo is free, is evolving rapidly and works reliably under most current operating systems.

A useful improvement to the language would be the addition of some form of structured variable, at least of arrays. For instance, the lack of arrays has lead, in the present project, to the choice of the second-order Runge-Kutta method over the better but more complicated fourth-order Runge-Kutta method.
Another nice improvement could be the implementation of some sort of a multi-windowed interface, so that different graphic views of the agent behaviour could be observed either alternatively or in parallel during the simulation run.
Finally, the developers should consider the addition of a grid to align the interface widgets and the possibility of opening several files in the same session.

CREDITS AND REFERENCES

[Maggi03] S. Maggi, in preparation
[Langenberg66] D. Langenberg et al., "The Josephson Effects", Scientific American, vol. 214, pp. 30-39 (May 1966).
[Anderson70] P. W. Anderson, "How Josephson Discovered His Effect", Physics Today, vol. 23, pp. 23-29 (1970).
[McCumber68] D. E. McCumber, "Effect of ac Impedance on dc Voltage-Current Characteristics of Superconductor
Weak-Link Junctions", Journal Applied Physics, vol. 39, pp. 3113-18 (1968).
[Kautz93] R. L. Kautz, "Chaos in a Computer-Animated Pendulum", Am. J. Phys., vol. 61, pp. 407-15 (1993)