Theoretical study of I-V characteristics in a coupled long Josephson junctions based on magnesium diboride superconductor

In the present work, the current-voltage (I-V) characteristics in a coupled long Josephson junction based on magnesium diboride are studied by establishing a system of equations of phase differences of various inter- and intra-band channels starting from the microscopic Hamiltonian of the junction system and simplifying it through the phenomenological procedures such as action, partition function, Hubbard-Stratonovich transformation (bosonization), Grassmann integral, saddle-point approximation, Goldstone mode, phase dependent effective Lagrangian and, finally, Euler-Lagrange equation of motion. The system of equations are solved using finite difference approximation for which the solution of unperturbed sine-Gordon equation is taken as the initial condition. Neumann boundary condition is maintained at both the ends so that the fluxon is capable of reflecting from the end of the system. The phase dependent current is calculated for different tunnel voltage and averaged out over space and time. The current-voltage characteristics are almost linear at low voltage and non-linear at higher voltage which indicates that the more complicated physical phenomena at this situation may occur. At some region of the characteristics, there exist a negative resistance which means that the junction system can be used in specific electronic devices such as oscillators, switches, memories etc. The non-linearity is also sensitive to the layer as well as to the junction thicknesses. Non-linearity occurs for lower voltage and for higher junction and layer thicknesses.


Introduction
Superconductivity of magnesium diboride (MgB 2 ) was discovered in 2001 with transition temperature of about 39 K [1]. Since its discovery as a superconductor, it has attracted the attention of the many researchers in the related fields because of its higher transition temperature than that of other metallic compounds. The two-gap nature of MgB 2 offers different types of physical phenomena which urges the researchers to work in the context of both theoretical and experimental prospects. The electronic structure of MgB 2 is similar to graphite which consists of honeycomb boron layers separated by magnesium layers [2]. The energy gaps are about Δ 1 = 2 meV which corresponds to two -bands and about Δ 2 = 7 meV which corresponds to -band. The state of Cooper pair corresponds to the gaps and are designated by the order parameters: 1 = Δ 1 e i 1 for the first gap and 2 = Δ 2 e i 2 for the second gap. The internal degree of freedom is the inter-band phase difference (ì, ) = 1 − 2 [3]. The degree of freedom can be increased by forming the stack of MgB 2 interlocked with the insulator such as SiO 2 , Al 2 O 3 etc., which is referred to as stacked Josephson junction. As a result, the Cooper pairs tunnel through the junction and the inter-band as well as intra-band phase textures are quite complicated.
In the present work, we derived a system of perturbed sine-Gordon equations for the coupled long Josephson junction. Starting from the microscopic BCS model Hamiltonian of the system and undertaking a number of steps of phenomenological path integral formalism, the phase dependent effective action and, hence, the effective Lagrangian density are derived. The system of equations of phase dynamics is derived by using Euler-Lagrange equation of motion to minimize the effective Langrangian density.
The system of equations are then solved numerically using the finite difference approximation imposing the Neumann boundary condition. The solution of unperturbed sine-Gordon equation is supplied as the initial condition. The computation is performed using OCTAVE 4.4 programming language. The present work is ended by giving the concluding remarks drawn from the computation.

Model Hamiltonian
The starting point of the present work is to write microscopic BCS Hamiltonian which is the total Hamiltonian of the system comprising the free Hamiltonian ( free ), pairing Hamiltonian ( pair ) and tunneling Hamiltonian ( T ) [4][5][6], that is Here, C † , , (ì, ) (C , , (ì, )) is the creation(annihilation) operator for fermion with spin = (↑ or ↓) for a given layer index and band index . These operators are the function of spatial coordinate ì and the imaginary time = −i . C † , , (ì, ) creates a fermion with spin at the given site (ì, ) and (ì, ) destroy the fermion therefrom. † (ì, ) and (ì, ) have the dimension of inverse square root of volume (i.e. Ω −1/2 ), with Ω being the total volume of the system. ì and 0 are the magnetic vector potential and electric scalar potential, respectively. * = 2 and is the electronic charge and is the mass of a fermion. The operator − ℏ∇ − * ì is called the canonical momentum operator.
Short-range or long-range phonon mediated fermions of oppsite spins form a syster of the pair of fermions. However, after the paring of such fermions, the fermionic nature is destroyed and a bosonic particle is formed. Such phonon mediated fermions having a bosonic property are called Cooper pair. , , is the coupling constant with the dimension of energy-volume (Jm 3 ). For the two-gap superconductor having s-and d-bands or is equal to ( , ). = refers to intra-band pairing and ≠ refers to inter-band pairing. Similarly, = refers to intra-layer and ≠ refers to inter-layer pairing.
The first term of equation (2.4) infers that a fermion of spin is destroyed in ( + 1) th layer and th band and is created in th layer and th band, while the second term is the complex conjugate of the first term. , +1 is the tunnel matrix element with the dimension of energy.

Action functional
According to the path-integral formalism of quantum mechanics, the action functional is defined as with L being the Lagrangian density. In terms of the total Hamiltonian, the action is defined as [5] Here, is the chemical potential, and is the total particle number, = 1 B , where B is the Boltzmann constant and is absolute temperature. is given as Now, the quantum mechanical partition function of the system can be written as Here, C is a column vector with elements C , , (ì, ), and C † is a row vector with elements C † , , (ì, ) while ∫ D [C † , C] represents the product of all integrals over the elements of C † and C.

Hubbard-Stratonovich transformation
The action functional associated with the pair Hamiltonian is in quartic form of four fermionic fields. Since ì and 0 are invariant under gauge transformation, the partition function of (2.9) can be rewritten as follows: Under the application of Hubbard-Stratonovich transformation, the quartic term of the pairing interaction can be reduced to quadratic, and the partition function becomes where,Δ(Δ) is the new fields which are bosonic in nature.Δ is a row vector containing the elements Δ , (ì, ) and Δ is a column vector containing the elements Δ , (ì, ). Applying special techniques of path integral formalism in this partition function and performing some matrix manipulation we could obtain the Lagrangian density as follows: At low temperature, the chemical potential is equal to the Fermi energy, i.e., = and = 0 since ↑ = ↓ . We also have, The effective Lagrangian is given by where, is the concentration of electronic charge, k F is the Fermi wave vector, TF = √︄ 0 π 2 ℏ 2 2 k F is the Thomas-Fermi charge screening length and L = √︂ 0 2 2 is the London penetration depth, ì and ì are electric and magnetic fields at layer and band.

13101-4
Theoretical study of I-V characteristics

Application to the long Josephson junction
Consider the stack of long Josephson junction with the lenght along -direction and junction system along -direction. External magnetic fields are applied along the -direction, which introduces a homogeneous phase difference along the -direction. The system is assumed to be uniform along the -direction and the problem becomes two dimensional. The electric field is along -direction. Now, the Lagrangian density in two-dimensional system becomes, as follows: where is the thickness of the superconducting layer and is the thickness of the junction material. is the dielectric constant of the junction material. The inter-band Josephson coupling constant is and Josephson tunneling coupling constant The -component of the electric field in between th and ( + ) th layer is and the -component of the magnetic field in between th and ( + ) th layer is We can introduce the gauge invariant phase difference , +1 as Then, we can have cos +1 − = cos , +1 and − = is the intra-layer inter-band phase difference. Hence, the equation (2.14) , +1 2

Coupled long Josephson junction system
In the coupled long Josephson junction, as shown in figure 1, there are eight channels for Cooper pair tunneling. The equations of phase dynamics can be obtained from the generalized equation , and so on.
In order to study the plasmon mode, the equation (2.23) can be linealized as   where¯and¯are the normalized frequency and wave vector, respectively.

Numerical computation and analysis
In order to perform the numerical computation, the equation (2.23) is discretized using the finite difference approximation. For this purpose, a uniform mesh in space and time is introduced with spacing and , respectively. At each point (¯,¯)

13101-8
Theoretical study of I-V characteristics where and are the total number of the points in space and time, respectively. The sine-Gordon equation is approximated by the second-order finite differences as where is the numerical approximation of the exact solution at (¯,¯). Applying this approximation, the perturbed sine-Gordon equation reads Providing the initial conditions to the junction system means supplying the initial information to the system at the starting time. In the present problem, the initial information is the kink (or anti-kink) solution of unperturbed sine-Gordon equation which can be generated by the appropriate electronic device which produces it as the trigger signal [7]. The solution of unperturbed sine-Gordon equation is where is the normalized speed of the kink ( = +1) or anti-kink ( = −1) and¯0 is its initial position. Hence, the initial condition for all channels of the junction system is The initial condition is approximated as (3.10) There are different boundary conditions that can be imposed on the system in order to control the state of kink or anti-kink. When (¯,¯) = 0 for¯= ± , then this condition will mirror the kink (anti-kink) and is known to be homogeneous Dirichlet boundary condition. The effect of moving will be demonstrated at the boundary¯= − to feed the domain with the incoming kink(anti-kink). The boundary condition then reads If the kink/anti-kink is to let reflecting from the boundary, then Neumann boundary condition,¯= 0 for¯= ± can be used [8].
In the present context, Neumann boundary condition is imposed which is approximated by central finite difference and yields 1 = −1 , at¯= − , and +1 = −1 , at¯= + . The Courant-Friedrichs-Lewy stability criteria, 2 2 < 1 should be maintained in order to get the stability of the kink/anti-kink solution. This criteria suggest us to take a very small time step as compared to the position step as far as possible. It is mandatory to pay the computational cost for the time steps to obtain the approximately calculated values, and reach a close agreement with those of theoretical values [9].
The most straightforward way to proceed the computational task is to introduce one array to hold at all¯at the time¯, a second array to hold all the at the time¯− 1 and a third array to hold a newly computed result at¯+ 1 . Then, it is looped through the code incrementing the time and shuffling the arrays appropriately [10].
The current can be calculated using [11][12][13] The current is averaged out over space and time as well as channel at different tunnel voltage which includes the element of equation (2.16) in the tunneling matrix. For the particular junction geometry, the tunneling matrix element is proportional to the bias voltage or tunnel voltage [14]. When the bias voltage is changed, then the tunnel matrix element also changes resulting in the change of the tunneling coupling constant of equation (2.16). This variation in the tunneling coupling constant significantly contributes to the soliton motion represented by the phase differences in various channels. The simulations were performed for a typical junction system of MgB 2 with superconducting layer thicknesses of 6Å, 9Å, 12Å and 15Å for different junction thicknesses 3Å, 6Å, 9Å and 12Å of SiO 2 . The junction and superconducting layer thicknesses are taken in the range of molecular dimension (i.e., Angstrom) as suggested by Giaever [15]. The dielectric constant of SiO 2 is taken as 3.7 and Fermi velocity of MgB 2 is taken as 4.7×10 5 m/s [16]. The computations were done using OCTAVE 4.4 programming language and the figures are generated using PYTHON3 in the Linux operating system on supercomputer platform.
To study the plasmon excitation, the dispersion relation defined by equation (2.25) is plotted and presented in figure 2. The figures 2a to 2d show that the band spectrum significantly depends on the junction and layer thicknesses. They appear in higher frequency states as the junction as well as the layer thicknesses are increased. The excited states can also be reached by increasing the bias-voltage that affects the tunneling matrix element. Hence, the tunneling coupling constant of equation (2.16) is also altered.
In order to study the soliton motion, the phase differences for all 8-channels are plotted against the normalized space and time. For this purpose, the LJJ of length 14 units was taken. The length is measured in the unit of inverse of Fermi wave vector (i.e. k −1 F ). The simulation was done up to the normalized time of 5 unit. Here, the time is measured in the inverse of Fermi frequency (i.e. −1 F ). From figures 3 to 6, it is observed that the initial kink in each channel greatly deforms its shape as the time lapses. As the kink reaches the boundary, it reflects back due to the application of Neumann boundary condition. As it reflects, there is a chance of producing the anti-kink or another kink. As a result, kink-kink or kink-anti-kink superposition may take place forming a complicated phase texture as time lapses. Some channels also show the collective behavior. The figures also show that the phase texture is highly sensitive to the junction and layer thicknesses as well as to the tunnel voltages. They become more complicated as the junction parameters are increased.
The I-V characteristics were studied by applying the voltage across the junction system. Since the tunneling matrix element is directly proportional to the bias voltage for the given junction and layer thicknesses [14], the tunneling coupling constant and hence the phase differences for all channels were computed for different applied voltages. Using the equation (3.13), the current for each channel can be computed at the given applied voltage as the function of space and time. The current is then averaged out over space and time as well as the channels for the given applied voltage and junction geometry. In this way, the current is calculated for each applied voltage ranging from 0 V to 1 V with the step of 0.01 V. The current is plotted against voltage as shown in figures 7a to 7d. The I-V characteristics for the junction thickness of = 3 Å are presented in figure 7a. The figure contains four curves for layer thicknesses of 6, 9, 12, and 15 Å. These curves seem to of be of non-ohmic nature with the existence of differential resistance. This resistive nature indicates that the applied voltage was consumed in order to proceed the tunneling of Cooper pairs. As the Cooper pairs reach the junction, they break into the normal charge carrier and they collide with the lattice in the junction. The junction system behaves as a conventional resistor. The graph shows that the positive differential resistance increases as the increment of the layer thickness. As the layer thickness increases, the population of Cooper pair also increases. For this reason, the collision frequency of the carriers was increased and a greater resistance persisted for the same applied voltage.
When the junction thickness was changed to 6 Å, an unusual type of I-V characteristics is obtained. For the given superconducting layer thickness, the I-V curve shows a positive differential resistance up to a certain applied voltage and then it shows a negative differential resistance. This peculiar behavior of the junction system indicates that there exists a complicated phenomenon. Due to the existence of a negative resistance at a certain voltage range, the device can be used as the energy stroage and oscillator. The origin of the noisy I-V characteristics may be interpreted in the following way. When the voltage is applied to the junction system, there are a lot of ways of dividing the value of this voltage into the voltages on those junctions. For this reason, there exist a lot of meta-stable states with different voltage distributions. Furthermore, it is quite possible that some meta-stable states are energetically very close to 13101-15 the stable state [11]. In such a case, the voltage distribution will be greatly changed and will cause a rapid oscillation of the dc current. For the junction thicknesses of 9 Å and 12 Å, the I-V curve showed even a complicated non-linear nature with the existence of a series of N-shaped differential resistances. The N-shaped differential resistance is the charanteristics of Gunn diode. Hence, the present junction system is also applicable to the low temperature electronic devices demanding the Gunn diode charactersitics. Another reason for the negative differential resistance is electromagnetic radiation due to the fluxonantifluxon (kink-antikin) transition between the plasmon excitation states during the soliton motion. As shown in the band spectrum (depicted in figure 2), the plasma wave is at the excited state causing the plasmon radiation. Therefore, the junction system in a particular voltage range can be used as a radiation chamber.

Conclusion
We conclude that the collision of fluxon and anti-fluxon as well as the in-phase or the out-phase of collective motion is more active for higher tunnel voltage. The current voltage characteristics are almost linear up to a certain tunnel voltage and then become non-linear. The non-linearity starts at a lower tunnel voltage for higher junction thicknesses as well as layer thicknesses. The linear region indicates that the junction system demonstrates a resistive nature while the non-linear condition confirms the existence of other complicated phenomena. Some nonlinear regions confirm the existence of a negative differential resistance so that the the junction system can be used in the electric devices that demand a negative resistance. One of the phenomena is the emission of electromagnetic radiation (e.g., microwave, THz etc.) due to the formation of meta-stable states as predicted by Koyama [13] and fluxon-antifluxon transition between them. The device can be used as a switching device, a memory device that operates in non-linear region. This might be the main region of THz radiation.