Exploration of fractional effects and a non-sinusoidal current-phase relation in Josephson junction arrays

Josephson junctions have a wide variety of applications, such as microwave oscillators, mixers and travelling wave power amplifiers (Kogan, 2022), (Guarcello et al., 2024) and a comb generator (Babenko et al., 2020). They are also of significant interest for quantum computations (Clarke and Wilhelm, 2008). Their dynamic behaviour has been a significant focus of research (e.g. Guarcello et al., 2023) and several recent research papers have explored the effect of noise on dynamical and transient behaviour (e.g. De Santis et al., 2022 and Fedorov and Pankratov, 2007). The paper considers a Josephson junction parallel array (Chevriaux et al., 2006) and explores the effect of a non-sinusoidal current-phase relation and spatial fractional derivatives on its dynamical behaviour. The closely connected sine-Gordon equation shall be considered to gain further insight into the behaviour. This equation models the behaviour of long Josephson junctions and Josephson junction arrays in the continuous limit.

The behaviour of Josephson junctions is described by partial differential equations involving the phase difference of the wave function across the junction between two superconductors. Of particular importance is the current-phase relation. Experimental techniques for its measurement are given in, for example, (Wetzstein et al., 2011) and (Golubov et al., 2004). It is generally assumed that Josephson junctions have a sinusoidal current-phase relation. However, as indicated in several works, for example, (Guarcello et al., 2025) and (Guarcello et al., 2024) and (Golubov et al., 2004), this is not necessarily the case and hence, the effect of a non-sinusoidal current-phase relation merits investigation. High-temperature superconductor and Fe-based superconductor junctions are examples of cases giving rise to an anharmonic current-phase relation (Askerzade and Salati, 2022). Graphene Josephson junctions have also been shown to have a highly skewed relationship (English et al., 2016). The authors in Canturk and Askerzade (2011) examined the influence of second harmonics on Josephson junction qubits and showed that phase qubits are strongly affected. The authors in Atanasova et al. (2011) examine the stability of the solutions of the double sine-Gordon model of the long Josephson junction. In Guarcello et al. (2025) the inclusion of a second-order harmonic affects the gain of a Josephson travelling-wave parametric amplifier and its stability. In Guarcello et al. (2024), a non-sinusoidal transparency-dependent current-phase relation is considered and the variation of the transparency is shown to affect the gain of the travelling-wave amplifier and the stability and route to chaos of the device.

In addition to this, the effect of fractional spatial derivatives shall be studied. Fractional spatial derivatives enable the modelling of nonlocal behaviour (Lischke et al., 2020). Fractional derivatives have been explored for several areas of science and electromagnetics, for example, (Antonini et al., 2021) and quantum mechanics (Bayin, 2016) and include global effects in system models. The accuracy of fractional time derivatives for modelling CMOS transmission line systems was shown experimentally in (Shang et al., 2013). Fractional effects have been examined in several previous works relating to the sine-Gordon model and Josephson junctions, for example (Bountis et al., 2025; Mares-Rincón et al., 2025; Ali et al., 2021; Macías-Díaz, 2017a; and Macías-Díaz, 2017b). There are various definitions of fractional derivatives (de Oliveira and Machado, 2014). In this work, the Riesz spatial derivative is selected as it reduces to the standard second-order finite difference for the second-order derivative. Its effect with a non-sinusoidal current-phase relation is considered.

The paper is concerned with the Josephson junction parallel array model (Chevriaux et al., 2006). This model becomes, in the continuous limit, the sine-Gordon equation when damping is ignored. The paper shall initially examine the effect of spatial fractional derivatives on the single sine-Gordon model, double sine-Gordon model and triple sine-Gordon model and compare the results with those obtained when the non-sinusoidal current-phase relation is used in the sine-Gordon model. The paper then proceeds to examine the Josephson junction parallel array model with integer and fractional spatial derivatives.

Results with an initial spatial distribution and without an input are compared for the various models. The application of the Josephson array as a sensitive signal detector is also examined. The non-sinusoidal current-phase relation and the spatial fractional derivatives are both seen to have a significant effect on the behaviour, with each changing the nature of the response.

The paper is concerned with a Josephson junction parallel array of $N$ junctions as described in (Chevriaux et al., 2006) and shown in Figure 1.

Let $JN={2,…,N-1}$

$un(t)$ is the phase difference across the junction $n$⁠. $Iin$is the input current to the array. $Iout$ is the output current through an output resistor, $Rout$⁠. $ς$ is a damping parameter.$k2$ is related to the Josephson inductance, $LJ$⁠. $Ls$ is the inductance of the wires connecting the junctions.

In the continuous limit, the model for the chain leads to the damped sine-Gordon Equation:

For integer-order spatial derivatives, when damping is ignored, exact solutions of the undamped sine-Gordon equation exist as detailed in many sources for example (Liu et al., 2006) and (Zun-Tao et al., 2005) and (Sun, 2015).

When $ς=0, k=1$ the travelling wave solutions take the form as follows:

These are known as kinks for the positive sign and anti-kinks for the negative sign. $v$ is the wave speed.

Another solution determined based on the work in (Sun, 2015) is as follows:

This solution is of interest as it is of similar form to a solution of the double sine-Gordon equation, which is of relevance when considering the inclusion of the second harmonic of the current-phase relation when modelling Josephson junctions.

The breather solution (Ferreira et al., 2008) takes the form of a soliton localised in space and oscillating in time as follows:

Its occurrence in long Josephson junctions is of significant research interest for potential applications in information processing (Fujii et al., 2008 and De Santis et al., 2022) and has been examined from a thermal transport perspective in De Santis et al. (2024). However, the experimental detection and stabilisation of breathers presents several challenges (e.g. De Santis et al., 2025).

Several past works, (e.g. Guarcello et al., 2025 and Canturk and Askerzade, 2011) have examined the effect of a second-order harmonic in the current-phase relation. In light of this, it is appropriate to consider the double sine-Gordon equation. This may be expressed as follows:

$A, β$ are constants.

Its exact solution may be determined as in (Sun, 2015):

$Q, k$ and $c$ are constants.

Non-sinusoidal current-phase relations for the Josephson junction of the form as follows:

or that scaled by its maximum value as follows:

have been considered in (Guarcello et al., 2024) which follows from that considered in (Golubov et al., 2004). $τ∈[0,1]$ is the transparency of the junction. $Ic$ is the critical current of the Josephson junction.

When $τ→0,$ the relationship approaches sinusoidal behaviour. When $τ→1$⁠, the relationship deviates significantly from a purely sinusoidal nature. Figure 2 compares the scaled function with $τ=0.1$ and $τ=0.9$⁠.

To assess the effect of the non-sinusoidal relation, consider a Fourier expansion of the non-sinusoidal relation:

As $τ→1$⁠, the magnitude of the higher harmonics increases. Figure 3 compares the magnitude of the harmonics for $τ=0.1$ and for $τ=0.9.$ In light of this, the sine-Gordon equation with multiple harmonics should be considered as follows:

For example, the triple sine-Gordon equation is as follows:

Fractional spatial derivatives have been in use in quantum mechanics, for example, in the fractional Schrodinger equation (Saha Ray, 2018). The fractional derivative enables inclusion of non-local effects and is suitable for modelling long-range spatial interactions (Kirkpatrick et al., 2013). Various numerical methods have been proposed, for example (Shen and Wang, 2024), (Owolabi and Atangana, 2016) and (Papadopoulos and Olver, 2024).

As in (Macías-Díaz, 2017b), the Riesz fractional derivative is considered for the present work. The single undamped sine-Gordon equation becomes as follows:

The Riesz fractional derivative may be expressed as follows for $0<α≤2, α ≠ 1$⁠:

$Γ$ is the gamma function.

For $α=1,$

$H$ is the Hilbert transform and p.v. is the principal value, (Macías-Díaz, 2017a).

In the Fourier domain, the Riesz fractional derivative is defined as follows:

$Fω$ is the Fourier transform of $ux,t$⁠. This property is used when performing the numerical simulations. The definition in (17) is valid for $0<α≤2$ (Bayin, 2016).

For the Josephson junction array, the discrete Riesz spatial fractional operator (Macías-Díaz, 2017b) is defined as follows:

$αj$ is the generalised binomial coefficient (Macías-Díaz, 2017a,b). When $α=2$⁠, the integer second-order central difference approximation of the spatial derivative is obtained. The discrete Riesz spatial operator with $0<α<2$ extends the second-order central difference to include global effects.

With the Riesz derivative, equation (1) becomes:

The solution of the sine-Gordon equation with multiple harmonics or with a different nonlinear current-phase relation may be determined numerically. In this contribution, the order-two Strang splitting method shall be used for numerical implementations. The splitting approach for solving sine-Gordon equations has been used in previous work, for example (Ham et al., 2022). A Fourier spectral method has been applied in (Zhen et al., 2026).

With Strang splitting, equations are separated into parts for which it is easier or more efficient to determine a solution. In relation to the sine-Gordon equations, it is appropriate to implement the spatial derivative in the Fourier domain and the nonlinear terms in the time domain using a numerical integrator such as the Runge-Kutta method. However, use of a Fourier basis does assume periodic solutions and this leads to oscillations in the solutions near the boundaries if boundary conditions are not periodic. Hence, the equations need to be modified to enable practical implementation with aperiodic boundary conditions.

Returning to equation (11), it is rearranged as follows:

and each part is solved separately. The first part is solved in the Fourier domain and the remaining part is solved in the time domain using the fourth-order Runge-Kutta method.

For practical implementations, the sine-Gordon equation is solved over a finite domain $[-L/2,L/2]$⁠. To enable use of the Fast Fourier Transform, zero boundary conditions should be considered. To enforce zero boundary conditions, the following modification is implemented:

$uL-vl-$ and $uL+vl+$ are the left and right boundary conditions.

Equation (20) is then transformed to as follows:

The Strang Fourier splitting method is applied to this equation and once $u∼(x,t)$ is determined, $ux,t$ follows from equation (21).

Consider the single sine-Gordon equation with an initial condition is as follows:

While theoretically, an infinite spatial domain is considered, for numerical computations $L$ is taken as 80 and results are shown for $x∈[-20,20]$⁠. The parameters are set as $k=1$ and $c=910$⁠. The timestep for simulations is set at 0.005.

Figure 4 shows the phase at $t=3$⁠, $u(x,3)$⁠, when the fractional spatial derivative $α$ is varied from 1.6–2. Note that $α=2$ corresponds to the standard second-order derivative and the exact solution for $α=2$ is also shown.

With the introduction of the fractional spatial derivative, there is a change in the rise time and the introduction of the oscillation is observed.

Now consider the case when a second harmonic is present. Consider equation (10) with $b1$ and $b2$ set equal to the first two Fourier coefficients of the Fourier series for the non-sinusoidal current phase in equation (9) with $τ=0.9$⁠. This value of $τ$ is selected to see the effect of a highly skewed current-phase relation. The initial condition remains the same.

Figure 5 shows the phase at $t=3,u(x,3)$⁠, when the fractional spatial derivative $α$ is varied from 1.6–2. The result with for the single sine-Gordon equation with $α=2$ is also shown. An oscillation is evident in all cases with the second harmonic including when $α=2$⁠. Comparing Figures 4 and 5, note that the second harmonic leads to oscillations of slightly greater amplitude – the maximum for $α=1.6$ increases from 4.03 to 4.13.

Now consider the triple sine-Gordon equation. $b1, b2$ and $b3$ are the Fourier coefficients determined for the non-linear current expansion in equation (9) with $τ=0.9$⁠. Figure 6 shows the results for the same initial condition. The extra harmonic results in increased oscillations. The maximum is 4.17 and the minimum of the $α=1.6$ oscillation is 1.18. Figure 7 shows the impact of fractional spatial derivatives on the result when the non-sinusoidal current-phase relation in equation (9) is used. While the minimum of the oscillations with the relationship in equation (9) has the lower value of 1.01, it is evident that the inclusion of three harmonics results in very similar oscillatory behaviour.

Consider the Josephson junction array modelled using equation (19) with the fractional derivative implementation of equation (18). $N=40$ sections are used for comparative purposes with the sine-Gordon results in the previous section.

Figures 8 and 9 show the results when the sinusoidal and non-sinusoidal [equation (9)] current -phase relations with $τ=0.9$ are used. The general form of the results is similar in nature to those in Figures 4 and 7 for the sine-Gordon equation. However, an oscillation is present in all results in Figure 8 unlike the $α=2$ case for the sine-Gordon equation. The finite length of the Josephson junction array results in the variation at the end points.

One application of the Josephson array is as a sensitive signal detector (Chevriaux et al., 2006). To examine this application and the effect of additional harmonics and fractional spatial derivatives, consider the same array with a damping coefficient $ς=0.1$ and an input $Iint=Asin⁡Ωt+s(t)$ where $s(t)$ is the signal to be detected.

As an example, let $A=0.5, Ω = 0.9$ and $st=0.2cosh⁡(0.1t-800).$

$s(t)$ corresponds to a pulse at $t=800$⁠. Figure 10 shows the input current, $Iin(t)$⁠. Figure 11 shows the voltage at node N for a sinusoidal current-phase relation. The region where the pulse is detected is shown in more detail in Figure 12. While the steady-state is greater for $α=2,$ the fractional derivatives delay the response to the pulse and dampen it.

Figure 13 and 14 consider the same example. However, this time three Fourier harmonics of the non-sinusoidal current phase [equation (8)] are included. In this case, the steady-state response is much greater and the occurrence of the pulse is less visible. The input frequency $Ω = 0.9$ is not in the forbidden gap. Note also the greater amplitude of the response as $α$ decreases. However if $Ω = 0.6$⁠, the result is as in Figure 15 and 16 and the pulse is much more visible.

When the complete non-sinusoidal function is considered, the results are similar. These results clearly indicate that the current-phase relation has a significant effect on the response of a Josephson junction array. The forbidden gap changes. In addition, the fractional effects change the damping of the response.

The paper has explored the effects of fractional spatial derivatives and non-sinusoidal current-phase relations with regard to Josephson junction arrays. The given results indicate that both features affect the behaviour of the systems. The sine-Gordon equation is initially considered, as it models long Josephson junctions and Josephson arrays in the limiting case. The additional harmonics and fractional derivative result in non-monotonic and oscillatory behaviour when the evolution of an initial spatial distribution is considered. Oscillatory responses are also observed in the Josephson junction array. One application of the Josephson junction is as a sensitive signal detector. However, non-sinusoidal current-phase relations affect the range of input frequencies in the forbidden gap and the transient behaviour of a Josephson array and would need to be considered in the design of applications. Fractional derivatives affect the damping of the response.

Future work could include examination of different forms of non-sinusoidal functions, different types of fractional derivatives and the inclusion of noise.