Superconductive Circuits and the General-Purpose Electronic Simulator APLAC

The general-purpose circuit simulator analysis program for linear active circuits (APLAC) is not widely known within the superconductor circuit community, regardless of its built-in Josephson junction model and capability of modeling the superconductive phase transition with controlled sources. We review the use of APLAC to model, e.g., noisy dc SQUIDs and transition-edge sensors. Based on an APLAC simulation, we also comment on the preference of voltage bias over current bias in dc SQUIDs in applications where large dynamic range and high bandwidth are simultaneously required.


I. INTRODUCTION
B EHAVIOR of linear superconductive circuits can be assessed with standard analytic techniques of the electronic circuit theory when the second Kirchhoff's law is replaced by the statement that rather than voltages, time integrals of voltages (i.e., fluxes) around closed loops sum to zero, or more generally to integer multiples of the flux quantum Φ 0 . Numerical treatment is however often necessary when second-order effects, such as the kinetic inductance or nonlinear phenomenon such as the Josephson effect or the superconductive phase transition, are of interest. In the early days of Josephson junction (JJ) studies, a lot of effort was invested in task-specific computer code for Josephson-dynamical simulation, but some general-purpose electronics simulators were also devised. These include WR-SPICE [1], PSCAN [2], and JSIM [3]. More recently, interest in large-scale implementations of RSFQ logic has driven further development of the electronic design automation (EDA) tools for Josephson circuits, leading to such simulators as [4] and [5]. Over the past decades, the slow speed of interpreted EDA simulators has prevented their serious use in studying, e.g., noisy Josephson dynamics, but the rapid increase of easily accessible computation power has made their use feasible in many cases.
Regardless of its moderately wide use at least within the Finnish superconductor community, in the literature one design tool has been rarely mentioned: the Analysis Program for Linear Active Circuits (APLAC). In this article, we describe how APLAC has been used at VTT for superconducting circuits, in particular for studying noise behavior of dc SQUIDs and dynamics of transition edge sensors (TESs).

II. APLAC AS A SIMULATION TOOL
Originating from the Radio Laboratory and the Circuit Theory Laboratory of the Helsinki University of Technology (HUT, currently Aalto University), roots of the APLAC program are almost as old [6], [7] as those of the more famous SPICE circuit simulator. From the mid-1980s onwards, a significant amount of resources was invested into APLAC development [8], owing to the interest of the Nokia Corporation to utilize it as a microwave design tool. A commercial company APLAC Solutions Corporation was founded in 1998 for further development, first acquired in 2005 by AWR Corporation, further acquired by the Cadence Design Systems recently.
It is moderately easy to implement the JJ by using controlled sources, with one voltage node reserved to represent the quantum phase [2]. An example JJ circuit implemented for LT-SPICE IV is available in [9]. A similar construction with controlled sources is also possible in APLAC. In the late 1980s, Heikki Seppä of VTT and Tapani Ryhänen of HUT, however, collaborated with Timo Veijola to build the JJ as a standard library item [10] compiled into the APLAC core. The addition of new elements was alleviated by the 1988 rewrite of the APLAC code in the C language with fully object-oriented techniques [11]. A builtin JJ model brought in faster simulation times as well as easy schematic capture and compact representation in the captured diagrams.
The JJ as a library component has continued to be available up to the current APLAC version 8.70 as of the year 2020. A quirk, worth noting in the JJ implementation in APLAC, is that the value of the flux quantum is precisely 2.07 mV·ps rather than its correct value 2.0678 … mV·ps.

A. DC-SQUID Example
An example of a current biased dc SQUID is shown in Fig. 1(a), as copy-pasted from the APLAC schematic editor. The schematic and netlist files are available in [9]. The SQUID loop is formed by inductors L1 and L2, coupled to the input coil L3 and Flux noise levels were computed by Fourier transforming the time trace of the output voltage in the manner detailed in [9]. (c) Voltage-biased SQUID characteristics were obtained with an added bias inductor [15].
(dimensionless junction capacitance) [12] β C = 0.7, and the dimensionless loop inductance of β L = 2L SQ I C /Φ 0 = 1.0. The time-domain simulation with Δt = 2 ps starts with a 2 14 -step initialization phase, whereby the SQUID bias and the input coil currents are driven from zero to their nominal values and the dynamics is settled. This is followed by a 2 17 -step stationary state, during which data are acquired, e.g., for noise analysis. There is a two-pole low-pass filter to remove the Josephson oscillation from the signal and to provide just the average voltage when plotting the SQUID characteristics. If characteristics alone are desired, the long steady state is not needed. Examples of the generated waveforms are found in [9].

B. Simulating the Noisy DC-SQUID
The noise estimation for a dc SQUID with APLAC can be performed in many ways. It is possible to use time-domain simulation to obtain the mixing coefficients, akin to the approach [13], as has been performed in [14]. As APLAC also includes the Harmonic Balance (HB) analysis mode, in the voltage-biased case (where the Josephson frequency is a priori known) it is in principle possible to separate the evolution of the quantum phase into the linearly progressing and periodic parts: where U B is the bias voltage. An HB circuit can then be constructed for χ(t) for obtaining mixing coefficients via simulation. These HB circuits have tended to be unreliable in convergence, however.
Most straightforward is the Langevin approach in the time domain, where time-domain Gaussian voltage noise generators are added in series with the junction shunt resistors. In the early versions of APLAC, we inserted the generators by hand as elaborated in [9], but in the more recent versions, we have found the built-in noise models of APLAC to function quite reliably, and the explicit generators are no longer necessary. The APLAC-internal noise models are activated by the Prepare … NOISE statement in the schematic Fig. 1(a).
In Fig. 1(b), the results from the current-biased circuit of Fig. 1(a) are shown when the circuit is invoked at 13 bias current and 100 applied flux values, and the output voltage noise during the 2 17 -step steady state is estimated by a Fourier transform of the time trace. In the noise simulation, the 2 14 -step initialization stage is followed by another 2 14 -step stage where a sinusoidal 10 mΦ 0p-p flux excitation is applied and the SQUID gain dV /dΦ is determined. Noise sources are then activated for 2 17 -step output voltage record. The output voltage noise is then scaled by the dV /dΦ to obtain the flux noise, which is plotted as contours in Fig. 1(b). This APLAC simulation consisting of 1300 iterations with 163 840 time steps at each iteration takes roughly 75 min on a Core-i7 equipped laptop computer running Windows 10. Fig. 1(c) shows the noise contours of a similar circuit with voltage bias. An inductor has been added between the bias voltage source and the SQUID, as described in [15]. The SQUID output current is recorded here and scaled by the dI/dΦ to obtain flux noise. Both current-and voltage-biased simulations were performed with the noise-generating temperature of T = 4.2 K.

C. Remarks on the Noisy DC-SQUID Results
Although the circuit of Fig. 1(a) implements a particular set of component values, the results can be scaled into dimensionless units [16], whereby they should be applicable to any dc-SQUID with the dimensionless junction capacitance β C = 0.7, the dimensionless loop inductance β L = 1.0, and the dimensionless noise energy Γ = 2πk B T /I C Φ 0 = 0.018. We are particularly interested to compare the numerically determined prefactor γ in the formula ε = γk B T L SQ C J for the SQUID energy resolution ε in terms of the loop inductance L SQ and junction capacitance C J . The flux noise contours in the plots Fig. 1(b) and (c) of Φ N = 0.4, 0.6, 1.0 and 1.5 μΦ 0 /Hz 1/2 correspond to γ = 9.4, 21, 60, and 130, respectively.
The canonical work for finding the SQUID noise optimum [16]- [18] has taken place solely with current bias, and only a narrow flux range of low-noise operating points has been deemed necessary as the flux-locked loop (FLL) [19] has traditionally been assumed to keep the SQUID flux close to the sweet spot. For instance, de Waal et al. [17] find γ = 12 in a narrow flux range with bias current I B = 1.6 · I C where not all flux settings lead to a finite SQUID voltage. We verify comparable or better γ values in the current-biased case [see Fig. 1 [20] is compatible with our data at I B > ∼ 2I C . In addition to the intrinsic flux noise, also the dV /dΦ finds a wider range of high values when using voltage bias [9]. The gain dV /dΦ determines how voltage noise of the readout amplifier translates into the effective flux noise. In the voltage biased case, an equivalent gain is determined as dV /dΦ = R D · dI/dΦ, where R D = dV /dI is the SQUID dynamic resistance at the setpoint.

(b)] in narrow ranges of applied flux when the bias current
We note that modern wideband SQUID applications, such as detector multiplexing [21], require a wide operating flux range, as large loop gain cannot be achieved in FLL-like negative feedback schemes if large bandwidth is desired [19]. Voltage bias [see Fig. 1(c)] not only yields more linear characteristics, but also a wider flux range of low-noise setpoints. Matched bias [22] lies between voltage and current bias and yields behavior (not shown here) between the current and voltage biased cases. In the matched bias, the load line due to the bias source is chosen roughly equal to the dynamic resistance of the SQUID.
The advantage of using an APLAC-like EDA tool is that it is trivial to augment the circuit [see Fig. 1(a)], e.g., by the backaction noise measurement or by adding parasitic circuit elements. Due to the targeted main commercial use of the simulator, APLAC natively contains a large number of microwave-related elements, such as transmission lines, stubs, and couplers.

D. Beyond the Resistively and Capacitively Shunted Junction (RCSJ) Model
It is possible to simulate the behavior of JJs shunted in a more complex manner than the standard RCSJ model. This opens the possibility of modeling, e.g., the DROS [23], the un-SQUID [24], or non-RCSJ style RSFQ circuits. Fig. 2(a) demonstrates a model that generates characteristics typically observed for a physical thin-film JJ. Here, the nonlinear subgap resistance is modeled by a nonlinear voltage-controlled current source, utilizing the hyperbolic tangent function. The current-biased model shows the hysteretic [12] IV-response [see Fig. 1(b)] familiar to any experimentalist. The response contains the zero-voltage branch when I B < I C , transition to the finite-voltage state at I B = I C , and the approach toward the zero voltage along the subgap branch when the bias current I B is lowered.
In this particular simulation, the bias current is designed not to return to zero at the end of the sweep, for the sake of demonstrating the retrapping behavior of hysteretic JJs. The current remains at I B = 1.6 μA for the final 2 ns at t > 8 ns [see the time-to-voltage plot Fig. 2)c)]. Such I B causes the average voltage in the subgap branch V JJ ≈ 200 μV to remain above the retrap voltage, and the supercurrent keeps oscillating at ∼100 GHz. If bias is lowered to I B = 1.5 μA, the supercurrent oscillation frequency gets sufficiently below the plasma frequency of the JJ that retrapping occurs, and the JJ decays to the non-oscillating zero-voltage state within the final 2 ns of the simulation.
Although the circuit [see Fig. 2(a)] has converged fine with APLAC 7.61, we note that the statement Prepare TRANMODE 0 needs to be added in APLAC 8.70 for stable simulation.  [25]. When the toggle input (blue, solid, voltage referred to the left Y-axis) is driven with Gaussian-shaped SFQ pulses, the Q-(dashed) and Q-hat (dash-dot) outputs show alternating SFQ pulses (referred to the right Y-axis). JJs are resistively shunted to obtain β C = 0.7.

E. More Complex Josephson Circuits
As a demonstration, we show in Fig. 3 the APLAC simulation of the RSFQ T-flip flop taken from [25]. Simulating over 2500 time steps with Δt = 0.1 ps consumes 0.2 s of CPU time on the Core-i7 laptop. The netlist and schematic are available in [9]. The waveforms indicate the correct T-flip flop operation, with alternating SFQ pulses appearing in the two outputs.

IV. SIMULATING TESS
Voltage nodes of a circuit simulator can be used to represent other physical quantities and current branches to represent other physical flows, in a similar manner as the quantum phase is represented as a node voltage in the LT-SPICE implementation of the JJ in [9]. Fig. 4 shows an implementation of a TES, taken from [26] and updated for the current APLAC version 8.70. We have used the TES model primarily to study the TES-readout interaction, e.g., electrothermal stability related to the LC resonator settling time in frequency-domain multiplexed SQUID readouts [21] and feasibility of the pulsed-bias time-domain multiplexing [26].
To get the TES simulation started, it is necessary to model both temperature-driven and bias current-driven (magnetic) transition from the superconducting to the resistive state. All variables in the simulation start from zero, and the TES is driven to its setpoint through the transition during the first few milliseconds of simulation. The model of the transition used here is The simulated TES current (plain line) as the bias voltage is swept from zero to 4 μV. The current shows the characteristic negative dynamic resistance, as well as the electrothermal oscillation due to the series inductance in the bias circuit at low voltages. The TES internal temperature (with markers) is also plotted. phenomenological, rather than one describing the accurate internal functioning of the TES.
The voltage of the "Ttes" node represents the TES temperature, and the node can be driven by current to model the heat flow due to incoming radiation. Examples of a dc-biased and ac-biased TES reacting to X-ray photons are given in [9], where a more detailed description of the model is also available.

V. CONCLUSION
We have found, over the years, APLAC to be a very useful tool in the design of practical SQUIDs and in understanding TES dynamics. Its long-standing commercial support and wide range of library components make it easier to simulate hybrids of Josephson-, TES-and more traditional electronic circuits. The APLAC facilities for hierarchical design alleviate modeling complex circuitry.
As a demonstration of APLAC capabilities, we have reproduced the canonical noise optima of the dc SQUID, while extending the set of simulated bias and flux region to exhaust the assumed "reasonable" range of values. In addition, we have extended the flux noise estimation to the voltage-biased case.

ACKNOWLEDGMENT
The author would like to thank J. Luomahaara and D. Hazra for commenting on the article.