[Rp] Reproduction of Step width enhancement in a pulse-driven Josephson junction

Introduction

A Josephson junction is a quantum mechanical device composed of two superconducting electrodes separated by a weak link [1]. For currents lower than a critical value $I_C$, coupled electrons (Cooper pairs) can cross the weak link without a potential difference (dc Josephson effect). When the current is increased above $I_C$, single electrons originated by the breakup of Cooper pairs begin to traverse the weak link. The potential difference $V$ between the two superconducting films becomes $\neq 0$ and a state is reached where the junction behaves as a resistance.

In modern Josephson junctions the weak link is usually a thin insulating tunnel barrier (SIS junction) [2], a normal metal film (SNS junction) [3] or a physical nanoconstriction (ScS junction) [4], [5]. Josephson junctions have found wide usage in several research fields, for example as building blocks for RSFQ digital electronics or quantum computers [6], or as radiation detectors and very sensitive magnetometers (SQUIDs) [7], [8], [9]. But the most successful application of Josephson junctions is surely in voltage metrology. A microwave radiation of frequency $f$ can phase lock the junction oscillations, producing the so-called Shapiro-steps, i.e., current steps at the quantized voltages $V_n$,

$$V_n = n \frac{h}{2 e} f, \quad n = 1, 2, ...$$

where $h$ and $e$ are the Plank constant and electron charge, respectively. The ac Josephson effect is at the basis of the current quantum voltage standard.

Besides its practical applications, the Josephson junction is important from a physical point of view because it has been the first device showing a quantum mechanical effect on a macroscopic scale.

An important research topic at the beginning of the ’90s was related to finding ways to increase the amplitude of the current steps induced by the microwave radiation (rf-induced steps). In fact, the stability of the lock between the phase of the junction and the applied microwave radiation – and therefore its insensitivity to noise events which might switch the junction from one quantized voltage to another, a crucial problem for voltage standard applications – is strongly dependent on the amplitude of the steps [10].

To increase the amplitude of the current steps, a non-sinusoidal microwave radiation may be used. In 1990 Monaco showed that, in the limit of a voltage-biased Josephson junction, adding together two phased microwaves of frequency $f$ and $2 f$ produces rf-induced current steps whose amplitudes are larger than those observed with a sinusoidal radiation [11].

Experiments on the so-called “biharmonic drive” readily confirmed these conclusions, albeit with some limitations due to the fact that the junctions could not realistically be considered as voltage biased [12], [13].

Extending further the idea, Monaco showed that, still in the limit of voltage bias, if the microwave radiation is composed of a train of delta functions, the rf-induced current steps could become as large as the critical current $I_C$.

However, a voltage bias configuration does not properly model a real Josephson junction, which should usually be considered as current biased. Also, a pulse train composed of delta functions is only a theoretical approximation and cannot be reproduced in actual experiments. This led to the idea to investigate what happened to a current-biased Josephson junction irradiated by a more realistic pulsed microwave signal [14], [15]. The reproduction of this investigation is the object of the present work.

Computational context

A first attempt to solve this problem was made using an electronic analog simulator of a Josephson junction [16], that could compute the relation between the applied current and the voltage ($I-V$ characteristic) of a current-biased junction, in the framework of the Stewart-McCumber RSJ junction model [17], [18]. The analog simulator was fast and simple to use, and could produce in just a few minutes on a Hewlett Packard 7475A2 pen plotter beautiful plots of the $I-V$ characteristics of the junction as a function of the simulated microwave signal. ¹

However, even if the electronic simulation was extremely fast, the analysis of the results required to measure by hand the amplitude of the rf-induced current steps visible on each $I-V$ characteristic, a tedious and error-prone task.

I then decided to develop a Fortran program to solve numerically the nonlinear second-order differential equation that models the Josephson junction [17], [18]. The idea was to calculate the $I-V$ characteristics of the junction as a function of the amplitude of the microwave signal, $\alpha_\mathrm{rf}$, for a given set of parameters characterizing the junction and the microwave, considering the three different cases of standard sinusoidal drive, biharmonic drive and pulsed drive.

To ease comparison of the results, normalized units were used throughout the calculations. The normalized junction voltage was $\eta$ and the normalized current $\alpha_\mathrm{dc}$. The main parameters of the simulation were: hysteresis parameter $\beta$, microwave frequency $\Omega$, amplitude of the microwave signal $\alpha_\mathrm{rf}$, pulse width $\rho$, integration time $\tau$. For a given set of junction parameters, usually $100$ different $I-V$ characteristics for increasing or decreasing values of $\alpha_\mathrm{rf}$ were calculated.

The first versions of the Fortran program were compiled under DOS 6.22 with Microsoft Fortran 5.1 and run on what was then a state-of-the-art PC, probably a Compaq Deskpro 486 with a math coprocessor, shared among several users of the lab.

The limitations of a PC for such a task soon become evident. A new simulation started automatically each evening and took the entire night to complete. People still using the PC late in the evening often inadvertently stopped the background process or simply shut the machine down without checking if there was another job running. At the end, I could run a full simulation only every two or three days.

After a few weeks of these mostly unsuccessful attempts, a colleague of another research group proposed me to use four DEC workstations running ULTRIX for my own simulations.² The machines were heavily used by his group during working hours, but sat mostly idle overnight. If I could manage to finish my runs before the start of the new work day, I was allowed to use this idle time for my own simulations. The colleague gave me a quick crash course on Unix and I was ready to go.

Porting my Fortran program from Microsoft Fortran 5.1 to ULTRIX was a breeze, and I quickly learnt how to use FTP to transfer the input configuration files and the output data files containing the results of the simulations back and forth from the DEC workstations to my Compaq 386 notebook, that was also my desktop computer.

Now each day I had four different sets of data files coming from the DEC workstations. Three of them were obtained by irradiating the junction with a pulsed drive with decreasing values of $\rho$, while keeping constant $\beta$ and $\Omega$. The fourth simulation was made by irradiating the junction with a sinusoidal drive, while keeping everything else equal. This last simulation was used as a reference, to compare the results obtained with a standard sinusoidal radiation with those obtained with progressively shorter pulses.

Each night I changed the values of $\beta$ or $\Omega$, to study the effect of these parameters on the behavior of the junction.

Again, the real problem was how to analyze all this data. A manual analysis like that needed with the electronic simulator was out of consideration. I decided to try the recently released Microsoft Visual Basic 1.0 for Windows, writing another program that calculated the size of the rf-induced current steps visible on the $I-V$ characteristics of the nightly simulations, as a function of $\alpha_\mathrm{rf}$.

The results of months of calculations were summarized in a paper published in the Journal of Applied Physics [14].

Digging into code

I like organisation, and I try keep all my past projects on my main workstation. Thus, finding the original source and data files of this project was only a question of locating the directory where the project was stored. Problems started to arise when I looked at the different files. The whole project was scattered into several directories, each containing many files with widely different names and dates. At first, trying to find an order in that chaos seemed impossible.

Normally I would have found all information needed in many notebooks full of detailed handwritten notes. Unfortunately, a couple of years ago most of my work notebooks were damaged by a water leak in the basement, and could not be recovered. The only option left was to check the files one by one.

After a thorough inspection of the whole project I recalled that: (1) my first attempts with the numerical simulations tried to use the more accurate McDonald-Johnson junction model [19], I later switched to the simplified Stewart-McCumber RSJ model because it was much faster and efficient in calculating the junction behavior [17], [18]; (2) file names attempted to reflect what the programs actually did, at least within the limits of DOS 8+3 naming scheme: in the same directory I could have a file ending with a “” that provided a textual output and another file ending with a “g” that gave a graphical output, and they differed only for a couple of DEFINESs that controlled the conditional compilation of the proper sections of the source code.

This multiplication of files might seem senseless today, when comfortable graphical interfaces, ultra-fast text editors with support of regular expressions and version control tools allow to change a large set of files in just a few seconds, but at that time it was probably the fastest, albeit very inefficient, way of working with source code; (3) the header of all source files contained detailed notes about the type of program, the compiler, the type of output and the dates of first and last revision of the source file, considerably easing the analysis of the different versions of the Fortran programs (Fig. 1); (4) the initial versions of the Fortran programs were monolithic, a single source file contained the whole code, that consisted in about 1.000 lines of Fortran. Only at a later time, better computing practices taught me to divide the monolithic code into multiple source files, compiled and linked together with a Makefile.

Another invaluable tool to analyze the different versions of the source files was Meld, an open source application available for all major operating systems that can perform a two- and even a three-way comparison of files and directories. Using Meld I quickly realized that the two most interesting source files were mcphase.for and mcp-work.for, both located in the mccumber/ directory (Fig. [fig:meld-comparison]).

The first program, mcphase.for, simulated the behavior of the junction for a single value of $\alpha_\mathrm{rf}$ read from the input configuration file mc-iv.dat, and saved the $I - V$ characteristic of the junction and its phase portrait (i.e., the relation between the phase and its time derivative, the latter being proportional to the junction voltage $V$) in the output file mc-iv.out.

Multiple calculations with several different value of $\alpha_\mathrm{rf}$ were performed by using a DOS batch file that basically choose one by one the configuration files containing the desired values of $\alpha_\mathrm{rf}$, renamed them to iv.dat, run the compiled executable mcphase.exe and at the end renamed the file containing the results, mc-iv.out, using a consistent naming scheme. I don’t recall why I choose this approach, but it was clearly very inefficient, as it required to prepare each day a long series of configuration files that differed only by the value of $\alpha_\mathrm{rf}$, and to update accordingly the DOS batch file that controlled the night calculations (Fig. [fig:batch-file]).

The second program, mcp-work.for was an improved version that could cycle across a set of several values of $\alpha_\mathrm{rf}$, producing a different output file for each value of $\alpha_\mathrm{rf}$. To simplify the later automatic analysis, it left out the phase portrait.

Clearly this was the program ported to the DEC workstations. Unfortunately, I could not find the actual source file used on these machines, perhaps because I worked directly on the workstations and never thought to copy back these files to my PCs. But Fortran is a very stable language and making mcp-work.for work on a modern machine was very easy.

As for the Microsoft Visual Basic 1.0 code, I found only two versions of the programs and the differences between them were minimal. Since both programs gave exactly the same results, I decided to stick with the version that had a still working precompiled binary file.

Porting Microsoft Fortran to modern Unix

Porting mcp-work.for to the XXI century so that it could be compiled with the modern open source and multiplatform gfortran Fortran compiler was very easy, thanks to the stability of the language across different versions and platforms. Only a few minor tweaks to the source code were needed.

All work has been done on macOS, which is essentially BSD Unix with a more appealing graphical interface, but it can be easily repeated on any modern Unix-like operating system such as Linux, and probably even on Windows, with the support of either the Windows Subsystem for Linux (for Windows 10) or of Cygwin (for earlier versions of the operating system).

Preprocessor directives

For reasons that go beyond my understanding Microsoft Fortran 5.1 did not use standard preprocessor directives, such as those supported by cpp or fpp, [20] but used a a slightly different proprietary syntax (Fig. [fig:preprocessor]). To support cpp, all was needed was to comment out all the $DEFINE directives in the header section of mcp-work.for and to replace the Microsoft Fortran 5.1 DEFINE blocks with standard cpp blocks throughout the code (Fig. [fig:preprocessor]).

The right directives are chosen now at compile-time. For example, the following command ³

$ gfortran -cpp -Dtextout -Dsingle -o mcp-work mcp-work.for

runs the cpp preprocessor before the gfortran compiler, selecting only the sections of code that produce a textual output (-Dtextout) and simulate the junction behavior with the sinusoidal drive (-Dsingle).

Filenames

Compilers based on Fortran 77, such as Microsoft Fortran 5.1, did not support dynamic memory allocation at runtime and required programmers to use fixed-length arrays and strings. Strings were used rather sparingly in Fortran code, so that was not a big deal. With an exception. My code defined the basename of all files as the 50-byte long character variable filename, attaching a proper extension to the input configuration file that contained the simulation parameters and to the output data files with the results of the calculations.

Under DOS that was not a problem, as DOS truncates file names to only 8 characters plus 3 characters for the extension, and excess characters were simply ignored. But under Unix file names have no practical limitations,⁴ and having all these 50 character-long filenames, mostly composed by blank characters, was ugly and complicated file management, in particular when using the command line interface. The solution was simple, as Fortran now has the TRIM function, that removes all trailing blanks from a string. Whenever a file is opened for reading or for writing, TRIM() is applied on-the-fly to the filename variable,

 OPEN (UNIT = 10, FILE = TRIM(filename)//'.dat', STATUS = 'OLD')

thus removing all extra blanks from the name of the file.

Edit descriptors

Microsoft Fortran 5.1 used the backslash (\) edit descriptor to prevent the addition of a line break at the end of a WRITE instruction. Modern Fortran compilers do not support this non-standard edit descriptor and return an error. This problem is avoided by removing the backslash from all WRITE instructions that include it.

Date and time

Microsoft Fortran 5.1 had two separate intrinsic subroutines to return the current date and time. In particular, CALL GETDAT(iyr, imon, iday), saved the date in the two-byte integer variables iyr, imon and iday, while CALL GETTIM(ihr, imin, imin, i100th) did the same for the current time, saving the return values in the integer variables ihr, imin, imin and i100th. The meaning of each returned variable should be self-explanatory.

Modern Fortran supports the single subroutine CALL DATE_AND_TIME(DATE, TIME, ZONE, VALUES), where all arguments are optional and can be specified by their dummy names (i.e., how Fortran calls the keyword arguments of a function call). In particular, DATE, TIME and ZONE are character variables, while VALUES is a one-dimensional array of 8 integers, where VALUES(1:3) corresponds to the year, month and day of the month, VALUES(4) is the time difference (in minutes) with UTC, and VALUES(5:8) are the hour, minute, second and milliseconds, respectively.

To minimize changes to the original source code, the calls to the GETDAT and GETTIM subroutines, were translated to a single call to DATE_AND_TIME, assigning the elements of the returned array of VALUES to integer variables named as in the original code (Fig. [fig:date-time]).

Compilation with gfortran

As noted above, mcp-work.for calculates the $I - V$ characteristics of the simulated junction for several different values of $\alpha_\mathrm{rf}$, producing one output file for each $I - V$ curve. The section of the code that defined the names of the output files was quite convoluted,

      il=0
      DO alpha_rf=0.0, 50.0, 0.5
      ...
c ------ define output file name(s)
         il=il+1
         il2=INT(il/100)
         il1=INT(il/10)-il2*10
         il0=il-il1*10-il2*100
         filewrite='PU'//CHAR(il2+47)//CHAR(il1+47)//CHAR(il0+47)
      ...
      END DO        ! repeat alpha_rf DO cycle
    

and used the integer variable il to count the cycle number, while the three integers il2, il1 and il0 contained the hundreds, tens and units digits of il, respectively. The CHAR function converts these integers to the corresponding ASCII characters, where ASCII code 48 corresponds to the 0 symbol and ASCII code 57 corresponds to 9. The name of the output file defined in the variable filewrite was built by concatenating a trailing constant string (’PU’ in the example above) to the three ASCII characters, using the double forward slash (//) operator.

I have no idea why I decided to build the name the output files in such a complex way, while it would have been much simpler to use the value of $\alpha_\mathrm{rf}$. In any case, it did not work under gfortran and prevented proper compilation of the code.

After some inspection it was apparent that the number 47 added to integer variables in the three CHAR function calls was the source of the error, and that it should be replaced by 48. The line defining the filewrite variable thus becomes,

      ...
      DO alpha_rf=0.0, 50.0, 0.5
      ...
c ------ define output file name(s)
         ...
         filewrite='PU'//CHAR(il2+48)//CHAR(il1+48)//CHAR(il0+48)
      ...
      END DO        ! repeat alpha_rf DO cycle

With this change mcp-work.for could compile flawlessly under gfortran. What is still puzzling is how the original line could work in Microsoft Fortran 5.1.

Visual Basic code

No attempt was made to try to run the original Visual Basic 1.0 program on a modern computer. Visual Basic is a dead language and long since has been replaced by Visual Basic .NET, which shares only the name with its forefather.

From the beginning, the only viable option to run a Visual Basic 1.0 application today was to rebuild the original development environment based on DOS 6.22 and Windows 3.11 in an emulator. The other possible alternative, try to setup an ancient PC still capable to run DOS and Windows 3.11, albeit in principle interesting to ensure a replication of the original paper at the hardware level, would have posed more problems than it solved, adding little to the accuracy of the reproduction itself.

I preferred to use the Parallels Desktop emulator for macOS, but popular alternatives such as VMware Workstation for Windows or the VirtualBox open source multi-platform emulator should work equally well.

I created a new empty virtual machine with minimal hardware requirements and installed in sequence DOS 6.22, Windows 3.11 and Visual Basic 1.0 (Fig. 2). While I was at it, and although I had already reproduced the Fortran part of the project, I also decided to install Microsoft Fortran 5.1 for DOS, to try to recreate as much as possible the original work environment.

All software packages were downloaded from the WinWorld web site, an invaluable resource for recovering old software packages. Even after so many years, the legitimacy of installing proprietary software in an emulator, might be questionable. But at the time I had regular licenses for all the above mentioned software and I guess to be at least morally authorized to continue to use those packages. Unfortunately, this also means that it is not possible to share the image of the virtual machine used to run the Visual Basic program on the paper’s GitHub repository, since it contains proprietary software.

The packages were originally distributed on several floppy disks, which had to be swapped whenever the installer required a new disk. The installation of these software packages in an emulator is close to how it was done back then. The only difference is that today the floppy disks are replaced by virtual file images and swapping disks is not done mechanically but requires to select a menu option in the emulator.

At the end of the installation process, the Windows 3.11 appeared as in Fig. 2. The default $640 \times 480$ pixel screen resolution of Windows 3.11 was woefully meager by today’s standards but, as I used the virtual environment almost exclusively to run the Visual Basic program, I didn’t bother to install the video drivers that could increase the screen resolution to a more comfortable $800 \times 600$ or $1024 \times 768$ pixel resolution.

Using the emulator required to transfer the source and binary Visual Basic files to the emulated DOS/Windows system. It is surely possible to make the emulated Windows 3.11 communicate with the host operating system through the network. But I found much easier to use again a virtual floppy disk to transfer all needed files from macOS to Windows 3.11 (and viceversa).

A new empty virtual floppy disk data.img can be easily created with the command-line utility dd available on macOS or Linux

dd if=/dev/zero of=data.img bs=1440k count=1

After creation, the virtual floppy disk must be mounted in the emulator and formatted under DOS or Windows 3.11 in the original MS-DOS FAT file system.

This step completed the preparation of the development environment, now it was time to test how all this behaved.

Running the programs

As already noted, mcp-work.for calculates the $I - V$ characteristics for a range of values of $\alpha_\mathrm{rf}$, saving each curve in a separate output file. The lower and upper limits and the step size of $\alpha_\mathrm{rf}$ are hardcoded in the Fortran source code ad every change requires a recompilation of mcp-work.for (a minor hassle, as the compilation takes just a couple of seconds on a modern machine).

Also the names of the output data files are hard-coded in mcp-work.for in the filewrite variable and are conventionally composed by a two-letter prefix (“SI” for the single drive, “BI” for the biharmonic drive and “PU” for the pulsed drive) followed by a three-digit integer that represents the cycle number (Section 4.5).

To avoid cluttering the mccumber/ directory that contains the Fortran source files with the output data files produced by the simulations, I created a new directory in the main project folder, 2020runs/, where I copied the mc-iv.dat configuration file needed to start the simulation (Fig. 1).

Running mcp-work.for requires three steps: compile mcp-work.for with the proper directives, switch to the 2020runs/ directory and run the mcp-work executable from there. The whole process is summarized below for the sinusoidal drive

$ gfortran -cpp -Dtextout -Dsingle -o mcp-work mcp-work.for
$ cd ../2020runs/
$ ../mccumber/mcp-work

The only modification needed to perform calculations using the pulsed drive is to change the -Dsingle directive to -Dpulsed

$ gfortran -cpp -Dtextout -Dpulsed -o mcp-work mcp-work.for
$ cd ../2020runs/
$ ../mccumber/mcp-work

On a recent (but not state-of-the art) machine the whole simulation with $100~\alpha_\mathrm{rf}$ steps takes around $5$ minutes for the single drive and $7$ minutes for the pulsed drive, and most of the time is spent printing on the terminal the calculated $I - V$ characteristics for each value of $\alpha_\mathrm{rf}$. Such feedback was useful at the time of the original calculation, as every new calculated point of the $I - V$ characteristics appeared on the screen after several tens of seconds, now the results scroll on the screen at a speed that makes them almost illegible. However, to keep as faithful as possible to the original project, I decided to continue to print the data points on the computer screen.

The $I - V$ characteristics calculated with the sinusoidal drive are shown in Fig. 3 for different values of $\alpha_\mathrm{rf}$. Without microwave radiation ($\alpha_\mathrm{rf} = 0$), the simulation produces the well-known $I - V$ characteristic of an overdamped Josephson junction (Fig. 3a), while for non-zero values of $\alpha_\mathrm{rf}$ the staircase-like structure of the rf-induced current steps appears on the $I - V$ curves (Fig. 3b and Fig. 3c).

For a pulsed drive, the $I - V$ characteristic without microwave radiation is identical to that calculated with the sinusoidal signal (Fig. 4a), while for $\alpha_\mathrm{rf} > 0.0$ the current steps induced by the pulsed drive are fewer than with the sinusoidal drive and can be nearly as wide as the critical current (compare Figs. 4a and 4c).

At the end of each run, the output files should be transferred to the Windows 3.11 virtual machine to be processed by stepampl, the Visual Basic application described in Section 5. However Windows 3.11 does not understand Unix line terminators and cannot read the output files without a preliminary conversion. The conversion can be easily done in the macOS Terminal by issuing the following command,

$ for f in $(ls *.out); do sed -i .bak s/$/$'\r'/ $f ; done

that changes the line terminators of the .out output files from the Unix format containing only a line-feed (LF) to the carriage return followed by a line-feed (CR-LF) format used by all versions of Microsoft Windows.⁵. The -i switch allows in-place conversion of each file. The original files are kept adding a .bak extension.

The output files in the proper Windows compatible format can now be transferred onto the virtual floppy disk image. When the transfer is done, the floppy disk image is unmounted from the host operating system and mounted in the virtual machine, making it visible to Windows 3.11. The output files are copied to an empty directory of Windows 3.11 and the Visual Basic application is started, either by running the precompiled stepampl.exe executable or by opening the Visual Basic project and running the program from there (Fig. [fig:stepampl]), saving the results in another text file with extension .STP (for steps) that could be transferred back to the host operating system via the virtual floppy disk image.

I also briefly tried to run the original Fortran code using the Microsoft Fortran 5.1 installed in the emulator. Compilation was fine but the resulting DOS program was extremely slow, taking about $30 - 35$ seconds for each $\alpha_\mathrm{rf}$ cycle and about $60$ minutes in total for the sinusoidal drive, more than a tenfold increase with respect to the native macOS version compiled with gfortran. Even considering the overhead of the emulator, the difference is too large not to be attributed to the low quality of the binary code generated by the Microsoft Fortran 5.1 compiler.

Results

At the time of writing the original paper the whole process had to be repeated each night for a different set of input parameters and for a different kind of microwave signal (single, biharmonic or pulsed).

Each night I used three of the available DEC workstations to simulate the junction behavior with the pulsed drive, using three different values of the (normalized) width of the pulse signal, $\rho$, while keeping constant all the other parameters, such $\beta$ and $\Omega$. The only other difference in these simulation was the range of variation of $\alpha_\mathrm{rf}$, which depended on the value of $\rho$ (shorter pulses require a much larger intensity of the rf-signal to have the same effect on the junction). The fourth workstation simulated the junction behavior with the standard sinusoidal drive, using exactly the same set of junction parameters.

The original paper contained all the information needed to reproduce the results shown in the figures, without having to repeat the whole analysis from scratch. The parameters used in all runs were: hysteresis parameter $\beta = 0.01$, frequency of the sinusoidal or pulsed rf signal $\Omega = 0.45$, current bias between $\alpha_\mathrm{dc} = -5.0$ and $\alpha_\mathrm{dc} = 5.0$, integration time $\tau = 500$, time step $\Delta \tau = 0.01$.

The simulation with the sinusoidal drive was performed by varying the amplitude of the microwave signal $\alpha_\mathrm{rf}$ between $0.0$ and $5.0$, with a step $\Delta \alpha_\mathrm{rf} = 0.05$. The three simulations with the pulsed drive were done using: (1) pulse width $\rho = 0.250$, $\alpha_\mathrm{rf} = 0.0 - 10.0$, $\Delta \alpha_\mathrm{rf} = 0.1$; (2) pulse width $\rho = 0.125$, $\alpha_\mathrm{rf} = 0.0 - 20.0$, $\Delta \alpha_\mathrm{rf} = 0.2$; (3) pulse width $\rho = 0.050$, $\alpha_\mathrm{rf} = 0.0 - 50.0$, $\Delta \alpha_\mathrm{rf} = 0.5$.

The resulting output and summary files were saved in separate folders in the 2020runs/ directory, named SINGLE/, PULS0250/, PULS125/, PULS0050/ after their DOS counterparts.

To reproduce the first and second figure of Ref. [14] I made the only concession to modernity. Instead of trying to recreate them with the plotting program used originally, probably Origin 2.0, I decided to write a couple of small R scripts that could automate the task. The results are shown in Fig. 5 and Fig. 6 and, as expected, are identical to those reported in the first two figures of Ref. [14]. The large vertical steps on the rightmost curves of Fig. 5 era the first rf-induced current steps, that can be nearly as large as the critical current $I_C$ without rf bias visible in the first $I - V$ characteristic on the left.

Figure 3 of the original paper could also be easily reproduced by plotting the maxima of the curves of Fig. 6, i.e., $\Delta i_n$ vs. $\alpha_\mathrm{rf}$, for the four different cases considered here.

Similar considerations can be made for the reproduction of Figure 4 of the original paper, which considers a slightly hysteretic junction with $\beta = 1.0$. For simplicity, I have chosen to show instead the $\Delta i_n$ vs. $\alpha_\mathrm{rf}$ curves for the two most significant cases of microwave signal, i.e., the sinusoidal drive and the pulsed drive with $\rho = 0.050$ (Fig. 7), from which the plots of Figure 4 could be easily replicated.

Code availability

All code used for this replication is available in the project’s GitHub repository. However, the term "all" should be taken with a grain of salt. While the Fortran is truly available to everyone and can be used as-is by compiling it with gfortran or with any other compatible modern Fortran compiler, the Visual Basic 1.0 code poses a completely different set of problems.

First of all, it can be run only by rebuilding an exact replica of the original development environment, as noted in Section 5. A task that without proper documentation can require a lot of trial-and error, as the author of the present paper discovered during the course of this work, when he found that the original code could not be imported in any later versions of Visual Basic for Windows and even in Visual Basic 1.0 for MS-DOS, released by Microsoft in parallel with the initial Windows version.

Second, the source code of a Visual basic program is a mix of Basic source files (having the usual .bas extension) and of Forms objects that compose the graphical interface of the program (with extension .frm), stored in some binary format and tied together in a container, known as a Visual Basic project (extension .mak), as shown in Fig. [fig:vb-project]). Making things worse, a Visual Basic project cannot be exported, the idea of code sharing or reuse across different applications was almost unknown at the time, at least in the commercial world.

The obvious consequence is that the reproduction described here would fail to comply, at least in part, with one of the basic requirements of the Ten Years Challenge, i.e. the availability of all the source code used for the reproduction.

But even if Visual Basic code seems locked for eternity in its proprietary binary format, there is a simple, albeit partial, solution to this problem, printing to a file. Printing was a true necessity back then to inspect or debug code. Computer screens were small and editors were primitive, the best way to have an overall view of a program during development was to print it on paper. Visual Basic 1.0 is no exception and can easily print both the Basic source code and the layout of each form composing the graphical interface. Once realized that, it is easy to set up a PostScript Printer in Windows 3.11 and to print the Basic sources and the forms to two separate PostScript files on the Windows 3.11 virtual disk. Transfer of these files to macOS is done again using the virtual floppy disk image and, once in macOS, it is a questions of seconds to convert the PostScript files to the more popular PDF format using the Preview application available in every version of macOS. The finishing touch is to use an R script (already developed for another project) to extract the text from the PDF file containing the Basic code, saving it to a true text source file, so that it can be easily inspected by everyone interested in this project.

Discussion

The main problem with this reproduction was the accidental loss of all my handwritten notes about the project, that compelled me to recover all information needed from the source and data files. Keeping good and updated documentation about any research project is paramount but equally important is to store all documentation in a safe location, where it can be easily retrieved.

Today most documents are in electronic form, making them even more prone to data loss whenever a disaster strikes. Implementing a good and reliable backup strategy on different forms of storage media should be a mandatory requirement of any research project, and this strategy should always combine local backups with backups on secondary storage locations, physically well separated from the original data source.

Starting the project today from scratch I would make very different choices about the programming languages to use for the simulation and data analysis steps.

Fortran, despite its venerable age, is still an excellent language for scientific programming, but today I would surely prefer Python, because of its flexibility, ease of use, and availability of excellent numerical libraries, such as NumPy and SciPy.

The main problem with Python is related to the tumultuous development of the language itself and of the thousands of available modules, which can cause incompatibilities even with code developed a few years ago. This problem can be, at least temporarily, be solved by using virtual environments, but a better standardized solution is strongly needed.

Another problem is related to its nature of interpreted language. However, the presence of well-documented Fortran and C bindings, can enhance performance of time-critical sections of code.

Starting today I would also avoid using a new language, as was Visual Basic 1.0 at the time for a scientific project. It is true that in 1993-94 it would have been very hard to foresee the rapid demise of Microsoft’s Visual Basic, nevertheless using a programming language only when its main core is stable, runs on a wide array of operating systems and is accepted by a wide community of developers is surely a safer bet.

Not every software project can afford to be as stable as TeX, the scientific typesetting system invented by the prominent mathematician and computer scientist Donald E. Knuth, that has reached a state where “it is unwise to make further improvements to the system [..] which should give the same results 100 years from now that they produce today” [21]. But on the other hand, a development environment that changes too much and too often or that is subjected to the whims of a single software company creates more problems than it solves.

Conclusions

Going back to my old paper has been an exceptionally interesting and instructive experience and I thank the organizers of this challenge for the opportunity offered.

However this is not only a nostalgic attitude. The reproducibility crisis is a serious issue today [22], that undermines scientific credibility and impacts the public’s trust in science, paving the way to all sorts of fake and unscientific beliefs. Being able to go back and reproduce what has been done in the past could ease the retraction of published papers containing fabricated, falsified, or modified data or results and could contribute to simplify the identification of future frauds.

This of course requires to share and make freely available the original data and the tools used to analyze them. A few years ago this requirement was impossible to fulfill in practice. The digital world in which we live makes it almost inevitable.

References