Coverage for qpdk / models / junction.py: 100%

1"""Josephson junction Models.""" 

2 

3import warnings 

4from functools import partial 

5 

6import jax 

7import jax.numpy as jnp 

8import sax 

9from sax.models.rf import admittance 

10 

11from qpdk.models.constants import DEFAULT_FREQUENCY, Φ_0 

12 

13 

14def _warn_if_overbiased(ib: float, ic: float) -> None: 

15 """Host callback to warn if bias current exceeds critical current.""" 

16 if jnp.any(jnp.abs(ib) >= ic): 

17 warnings.warn( 

18 "DC bias |I_b| >= I_c detected. Linearized RCSJ model is invalid in the voltage state. " 

19 "Please check your 'ib' and 'ic' ('ic_tot', 'flux', 'asymmetry') inputs.", 

20 RuntimeWarning, 

21 stacklevel=0, 

22 ) 

23 

24 

25@partial(jax.jit, inline=True) 

26def josephson_junction( 

27 *, 

28 f: sax.FloatArrayLike = DEFAULT_FREQUENCY, 

29 ic: sax.Float = 1e-6, 

30 capacitance: sax.Float = 5e-15, 

31 resistance: sax.Float = 10e3, 

32 ib: sax.Float = 0.0, 

33) -> sax.SDict: 

34 r"""Josephson junction (RCSJ) small-signal Sax model. 

35 

36 Linearized RCSJ model consisting of a bias-dependent Josephson inductance 

37 in parallel with capacitance and resistance. 

38 

39 Valid in the superconducting (zero-voltage) state and for small AC signals. 

40 

41 Default capacitance taken from :cite:`shcherbakovaFabricationMeasurementsHybrid2015`. 

42 

43 See :cite:`McCumber1968` for details. 

44 

45 Args: 

46 f: Array of frequency points in Hz 

47 ic: Critical current :math:`I_c` in Amperes 

48 capacitance: Junction capacitance :math:`C` in Farads 

49 resistance: Shunt resistance :math:`R` in Ohms 

50 ib: DC bias current :math:`I_b` in Amperes (:math:`\|I_b\| < I_c`) 

51 

52 Returns: 

53 sax.SDict: S-parameters dictionary 

54 """ 

55 jax.debug.callback(_warn_if_overbiased, ib, ic) 

56 

57 ω = 2 * jnp.pi * jnp.asarray(f) 

58 

59 # Bias-dependent phase factor 

60 # jnp.clip ensures we don't get NaNs during compilation tracing or slightly overbiased states 

61 cos_Φ_0 = jnp.sqrt(jnp.clip(1.0 - (ib / ic) ** 2, a_min=1e-10)) 

62 

63 # Josephson inductance 

64 Lⱼ = Φ_0 / (2 * jnp.pi * ic * cos_Φ_0) 

65 

66 # Admittances (parallel RCSJ) 

67 Y_R = 1 / resistance 

68 Y_C = 1j * ω * capacitance 

69 Y_L = 1 / (1j * ω * Lⱼ) 

70 

71 # Total admittance 

72 Y_JJ = Y_R + Y_C + Y_L 

73 

74 return admittance(f=f, y=Y_JJ) 

75 

76 

77@partial(jax.jit, inline=True) 

78def squid_junction( 

79 *, 

80 f: sax.FloatArrayLike = DEFAULT_FREQUENCY, 

81 ic_tot: sax.Float = 2e-6, 

82 asymmetry: sax.Float = 0.0, 

83 capacitance: sax.Float = 10e-15, 

84 resistance: sax.Float = 5e3, 

85 ib: sax.Float = 0.0, 

86 flux: sax.Float = 0.0, 

87) -> sax.SDict: 

88 r"""DC SQUID small-signal Sax model in the zero-screening limit. 

89 

90 Treats the DC SQUID as a single effective RCSJ whose critical current is 

91 tunable by an external magnetic flux. Assumes negligible loop geometric inductance. 

92 

93 See :cite:`kochChargeinsensitiveQubitDesign2007a` and :cite:`tinkhamIntroductionSuperconductivity2015` 

94 for details on asymmetric SQUIDs and effective Josephson inductance. 

95 

96 Args: 

97 f: Array of frequency points in Hz 

98 ic_tot: Total critical current sum :math:`I_{c1} + I_{c2}` in Amperes 

99 asymmetry: Junction asymmetry :math:`(I_{c1} - I_{c2}) / I_{c,tot}` 

100 capacitance: Total SQUID capacitance :math:`C_1 + C_2` in Farads 

101 resistance: Total SQUID shunt resistance :math:`R_1 || R_2` in Ohms 

102 ib: DC bias current :math:`I_b` in Amperes 

103 flux: External magnetic flux :math:`\Phi_{ext}` in Webers 

104 

105 Returns: 

106 sax.SDict: S-parameters dictionary 

107 """ 

108 ω = 2 * jnp.pi * jnp.asarray(f) 

109 

110 # Normalized flux 

111 phi_reduced = jnp.pi * flux / Φ_0 

112 

113 # Flux-tunable critical current for an asymmetric SQUID 

114 ic_squid = ic_tot * jnp.sqrt( 

115 jnp.cos(phi_reduced) ** 2 + (asymmetry**2) * jnp.sin(phi_reduced) ** 2 

116 ) 

117 

118 jax.debug.callback(_warn_if_overbiased, ib, ic_squid) 

119 

120 # Bias-dependent phase factor 

121 # Guard against ic_squid == 0 (e.g. symmetric SQUID at half-flux quantum) 

122 ratio_sq = jnp.where( 

123 ic_squid > 0, (ib / jnp.where(ic_squid > 0, ic_squid, 1.0)) ** 2, 0.0 

124 ) 

125 cos_Φ_0_eff = jnp.sqrt(jnp.clip(1.0 - ratio_sq, a_min=1e-10)) 

126 

127 # Effective Josephson inductance (will go to inf when ic_squid == 0) 

128 Lⱼ_eff = Φ_0 / (2 * jnp.pi * jnp.where(ic_squid > 0, ic_squid, 1e-20) * cos_Φ_0_eff) 

129 

130 # Admittances (parallel RCSJ) 

131 Y_R = 1 / resistance 

132 Y_C = 1j * ω * capacitance 

133 # Guard against division by zero in Y_L if L_j_eff somehow goes to zero or Y_L -> 0 if L_j_eff is inf 

134 Y_L = jnp.where(ic_squid > 0, 1 / (1j * ω * Lⱼ_eff), 0j) 

135 

136 # Total SQUID admittance 

137 Y_SQUID = Y_R + Y_C + Y_L 

138 

139 return admittance(f=f, y=Y_SQUID) 

140 

141 

142if __name__ == "__main__": 

143 import matplotlib.pyplot as plt 

144 

145 sax.set_port_naming_strategy("optical") 

146 

147 f = jnp.linspace(1e9, 10e9, 500) 

148 s = josephson_junction(f=f, ic=1e-6, capacitance=50e-15, resistance=10e3, ib=0.5e-6) 

149 plt.figure() 

150 plt.plot(f / 1e9, jnp.abs(s[("o1", "o1")]), label="|S11|") 

151 plt.plot(f / 1e9, jnp.abs(s[("o1", "o2")]), label="|S12|") 

152 plt.plot(f / 1e9, jnp.abs(s[("o2", "o2")]), label="|S22|") 

153 plt.xlabel("Frequency [GHz]") 

154 plt.ylabel("Magnitude") 

155 plt.legend() 

156 plt.title("Josephson Junction S-parameters") 

157 plt.show()