Numerical study for electromagnetic wave emission in thin samples of intrinsic Josephson junctions

Physica C 471 (2011) 1202–1205

Contents lists available at ScienceDirect

Physica C journal homepage: www.elsevier.com/locate/physc

Numerical study for electromagnetic wave emission in thin samples of intrinsic Josephson junctions T. Koyama a,c,⇑, H. Matsumoto a,c, Y. Ohta b,c, M. Machida b,c a

Institute for Materials Research, Tohoku University, Sendai 980-8577, Japan CCSE, Japan Atomic Energy Research Agency, Tokyo 110-0015, Japan c JST-CREST, Kawaguchi 332-0012, Japan b

a r t i c l e

i n f o

Article history: Available online 14 May 2011 Keywords: Intrinsic Josephson junctions THz emission

a b s t r a c t Emission of THz electromagnetic waves from thin samples of intrinsic Josephson junctions (IJJ’s) is numerically studied, using the xz-model. We show that the spatial symmetry of the electromagnetic excitations corresponding to the p-cavity mode is different from that of the 2p-cavity mode in the IJJ’s where the junction parameters such as the Josephson critical current are weakly inhomogeneous. In such IJJ’s the emission in the [0 0 1] direction, which is forbidden in the dipole emission, appears at the p-cavity mode resonance, whereas it is not observed in the 2p-cavity mode resonance. It is also shown that the strong emission occurs when the transition between branches in the I–V characteristics takes place. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction

2. Coupled dynamical equations in the xz-model

Since the discovery of the THz wave emission from Bi-2212 intrinsic Josephson junctions (IJJ’s) [1] a lot of papers exploring the physics of electromagnetic wave (EM) emission from IJJ’s have been published [2–11]. In this paper we investigate numerically the EM wave emission from thin samples of the IJJ’s on the basis of the xz-mode, i.e., a 2D multi-junction model, taking account of the inductive and capacitive couplings between junctions. In our numerical simulations the EM field excited in the voltage state of the IJJ’S is calculated in the whole space composed of IJJ’s, electrodes and the vacuum outside the junctions simultaneously without using the dynamical boundary condition at the junction edges [12], which is used in many theoretical works. Then, one can obtain the power of emitted EM waves directly in the vacuum region in our numerical method. A detailed numerical study on the EM wave emission generated by the in-phase motion of the phasedifferences in IJJ’s, using the xz-model, is presented in our previous paper [10]. In this paper, we extend the calculations to more general cases in which the degree of freedom of the out-of-phase motion of the phase-differences is incorporated. We clarify the EM modes excited at the cavity-mode resonance and the power distribution of the emitted EM waves. It is also shown that the strong emission takes place in the retrapping region where the uniform branch in the I–V characteristics becomes unstable.

In this paper we investigate the electromagnetic (EM) wave emission from a thin sample of intrinsic Josephson junction stacks (IJJ’s) in the voltage state under an applied DC bias current. Consider the IJJ’s composed of N junctions covered with normal-metal electrodes as shown in Fig. 1. We suppose that the thickness of the IJJ’s, Dz = N(s + d) with s and d being the thicknesses of the superconducting and insulating layers, is much thinner than the wavelength of emitted EM waves, e.g., Dz  1 lm. In this paper the xz-model, which is a 2D model presented in Ref. [10], is used for investigating the electromagnetic excitations and the THz wave emission in the IJJ’s. This model describes the dynamics of the EM field in the whole system, i.e., the IJJ’s, electrodes and the vacuum. In the region inside the IJJ’s the gauge-invariant phase differences, /(x, ‘, t), and the electric and magnetic fields, Ex(x, ‘, t), Ez(x, ‘, t) and By(x, ‘, t), with ‘ being the junction index (‘ = 1, . . . , N) in the xz-model satisfy the coupled dynamical equations as follows:

 2ed  1  gD2z By ðx; ‘; tÞ; hc 2ed ð1  aD2z ÞEz ðx; ‘; tÞ; @ t /ðx; ‘; tÞ ¼ h

rx /ðx; ‘; tÞ ¼

0921-4534/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.physc.2011.05.159

ð2Þ

1 4p @ t Ez ðx; ‘; tÞ ¼ rx By ðx; ‘; tÞ  ðj sin /ðx; ‘; tÞ þ rEz ðx; ‘; tÞÞ; c c c ð3Þ

⇑ Corresponding author at: Institute for Materials Research, Tohoku University, Sendai 980-8577, Japan. Tel.: +81 22 215 2008; fax: +81 22 227 1469. E-mail address: [email protected] (T. Koyama).

ð1Þ

Ex ðx; ‘; tÞ ¼ 0;

ð4Þ

T. Koyama et al. / Physica C 471 (2011) 1202–1205

vacuum

I

electrode

IJJ's D=N(d+s)

I

electrode

Fig. 1. Schematic view of the xz-model. The IJJ’s covered with normal electrodes are located in the vacuum. A bias current I is injected into the electrodes as a current beam as shown by bold arrows.

with  being the permittivity of the insulating layers. In Eqs. (1) and (2), which are known as the generalized Josephson relations, g and a are the inductive and capacitive coupling constants, respectively, and D2z is the 2nd-rank difference operator, i.e., D2z f‘ ¼ f‘þ1  2f ‘ þ f‘1 . Eq. (3) is one of the Maxwell equations, in which the tunneling current across the junctions is given by the sum of the Josephson current with the critical current density jc and the quasi-particle tunneling current with conductivity r. Note that the electric field and the current in this region are assumed to have only the components perpendicular to the junction planes. In the following the electrodes attached to the IJJ’s are assumed to be made of a normal metal which satisfy Ohm’s law with conductivity rL in this paper. Then, the Maxwell equations inside the electrodes take the form as

1 4p @ t L Ex ðx; z; tÞ ¼ rz By ðx; z; tÞ  ðj ðx; z; tÞ þ rL Ez ðx; z; tÞÞ; c c ext

ð5Þ

1 @ t L Ex ðx; z; tÞ ¼ rz By ðx; z; tÞ; c

ð6Þ

1 @ t By ðx; z; tÞ ¼ rx Ez ðx; z; tÞ; c

ð7Þ

with L being the permittivity of the electrodes. In Eq. (5) jext(x, z, t) represents an injected current beam, which is regarded as an applied DC bias current. In the vacuum region outside the sample the Maxwell equations are given by

1 4p @ t 0 Ex ðx; z; tÞ ¼ rz By ðx; z; tÞ  j ðx; z; tÞ; c c ext 1 @ t L Ex ðx; z; tÞ ¼ rz By ðx; z; tÞ; c 1 @ t By ðx; z; tÞ ¼ rx Ez ðx; z; tÞ: c

ð8Þ ð9Þ ð10Þ

In Eqs. (5) and (8) the spatial distribution of jext(x, z, t) is chosen as

jext ðx; z; tÞ ¼ jext hðx þ Lx =2Þhðx  Lx =2Þ  ½hðW z =2  zÞÞ þ hðz  W z =2Þ;

ð11Þ

where Wz is the penetration depth of the current beam. We have checked that the numerical results are not sensitive to the value of Wz. Note that from Eqs. (1)–(3) it follows the well-known coupled sine–Gordon type equations,

h

1  aD2z 

¼

1

i1

c2

! @ 2t /ðx; ‘; tÞ

3. Numerical results for the emitted power and the I–V characteristics

1

k2J r2x /ðx; ‘; tÞ  sin /ðx; ‘; tÞ h i1 pffiffiffik  c @ t /ðx; ‘; tÞ:  b 1  aD2z c

g

 D2z

k2c

1203

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi pffiffiffi where kc ¼  hc2 =8pedjc ; kJ ¼ kc = g and b ¼ 4pr= ckc . The coupled Eq. (12) have two spatial scales that rule the spatial variation of the phase differences in the in-plane direction, i.e., the c-axis penetration depth kc and the Josephson length kJ. The lengths, kc and kJ, give the spatial scales of the in-phase and out-of-phase components of the excited EM field, respectively. In Bi-2212 IJJ’s we find kc P 100 lm and kJ  1 lm, namely, g P 104. Note that the two scales are greatly different in magnitude in Bi-2212 IJJ’s in which the THz wave emission takes place. The electromagnetic excitations of both spatial scales in the in-plane direction are possible inside the IJJ’s. However, the excitations of the scale kJ with frequencies in a THz range are confined inside the IJJ’s, that is, the electromagnetic excitations that cause THz wave emission in Bi-2212 IJJ’s are long wavelength in-phase modes with the scale kc, which has been confirmed in experiments [1,3]. The numerical method to solve the dynamical equations in the IJJ’s in the whole space, which is based on the ‘‘Finite Difference Time Domain method’’, was given in our previous paper [10]. In that paper we solved Eqs. (1)–(10), assuming synchronized motion for the gauge-invariant phase differences, i,e., ‘-independent /(x, ‘, t), and discussed the THz wave emission originating from the in-phase modes. In this paper we extend the calculations to the ‘-dependent general case in which the effect of multi-junctions arising from the inductive and capacitive couplings between junctions is explicitly taken into account. In the IJJ’s each junction has a nano-scale thickness, for example, s + d  15 Å in Bi-2212 IJJ’s. As a result, the EM field excited by non-synchronized motion of the gauge-invariant phase differences varies on a scale of nanometers along the c-direction inside the IJJ’s. However, the EM oscillations with a nano-scale wavelength cannot excite EM waves in a THz range that can reach the far-field region in the vacuum, which indicates that only the in-phase components included in the excited EM modes in the IJJ’s can be the source of EM waves outside the IJJ’s. It is also noted that the wavelength of the EM wave with a frequency in a THz range is on the order of 102 lm in the vacuum except for the near-field region, that is, the spatial scales of the EM field along the c-direction are greatly different between inside and outside the IJJ’s. Noting these facts, we perform numerical calculations in the two regions inside and outside the IJJ’s separately with different spatial scales along the c-direction and combine the EM fields in the two regions at the junction edges under the condition that the uniform component of the EM field inside the IJJ’s connects continuously to the EM field in the vacuum outside the IJJ’s. In this paper, we further use an approximation to make numerical calculations performable for the IJJ’s with a thickness much thinner than the wavelength of emitted EM waves, kEM, i.e., Dz  1 lm  kEM  102 lm, which corresponds to the experiments for Bi-2212 mesa samples. In our numerical calculations the size of the meshes, Dx  Dz, in the vacuum on which EM field is calculated is taken to be larger than the thickness of the IJJ’s, for example, Dx  Dz = 0.01kc  0.01kc and Dz < 0.01kc. This situation is realized when Dx = Dz = 0.01kc = 2 lm, Dz ’ 1 lm and kc = 200 lm. In this case the junction edges on the left and right sides are reduced to two points on the meshes. The values of the EM field on these points are taken to be equal to the mean values of the EM fields at the junction edges in the IJJ’s. The details of our numerical method used in this paper will be given in a separate paper.

ð12Þ

We solve Eqs. (1)–(3) for the IJJ’s composed of 20 junctions under the periodic boundary condition, h‘ + 20(x, t) = h‘(x, t). The in-plane size of the junctions is chosen as Lx = 0.5kc, as an example, in this paper, which leads to the cavity mode frequencies, xn/

T. Koyama et al. / Physica C 471 (2011) 1202–1205

pffiffiffi xp = 2np (n = 1, 2,. . .), with xp  c= kc being the Josephson plasma frequency, and the capacitive and inductive coupling constants are taken as a = 0.1 and g = 104, which are allowed for Bi-2212 IJJ’s. In the following numerical calculations we introduce tiny inhomogeneity of the Josephson critical current density along the in-plane direction, assuming the functional form, jc(x)/jc = 1  jc0[1  tan h ((x + Lx/2)/wc0)], where the values of the parameters, jc0 and wc0, are chosen as jc0 = 0.001 and wc0 = 0.05kc in this paper. One may expect that the system can easily reach the most stable steady state in the presence of inhomogeneity. Let us now present the numerical results. Fig. 2 shows the voltage dependence of the power (S = jE  Bj/4p) emitted in the [1 0 0] direction, which is measured at a distance of 1.6kc from the center of the IJJ’s in the vacuum region. As seen in this figure, we have two peaks in this voltage region, though the peak located at a higher voltage value is very weak. Note that the voltages at which the peaks appear are nearly V/Vp = 2p and 4p with Vp = ⁄xp/2e, respectively. Since the Josephson frequencies at these voltages are close to the cavity mode frequencies, x1 (p-mode) and x2(2p-mode), in the present IJJ’s, these two peaks are understood to originate from the cavity mode resonance. Note that the p-mode is the half-wavelength nodeless mode and the 2p-mode is the one-wavelength mode with a single node. In Fig. 3 we plot the power emitted in the three directions, [0 0 1], [1 0 1] and [1 0 0], at these two resonances. It is seen that the directions in which the strongest emission occurs are different between these two resonance modes. The analysis for the EM field excited inside the IJJ’s indicates that the spatial pattern in the ac EM field is symmetric, i.e., Ez(x, ‘, t) = Ez(x, ‘, t) and By(x, ‘, t) = By(x, ‘, t), at the resonance with the 2p-mode, which causes a dipole-like emission, i.e., S / sin2h with h being the angle from the c-direction. On the other hand, the emission at the p-mode resonance is not the dipole-like, that is, as seen in Fig. 3, the strongest emission takes place in a direction between [0 0 1] and [1 0 0] and also the emission in the [0 0 1] direction does not vanish. Our numerical result for the spatial pattern of the ac EM field excited at the p-mode resonance has both symmetric and anti-symmetric components. The emission pattern at the p-mode resonance given in Fig. 3 is consistent with the experiments [3]. Fig. 4 shows the I–V characteristics in the voltage region near the p-mode resonance. In this figure we also plot three lines satis-

(×10-5)

2.5 (×10-7)

(×10-5) 2.5

[101] 2

2

[001]

[100]

π -mode

1.5

1.5

S/Sp

1204

1

1

[100] [101]

0.5

0.5

2 π-mode [001]

0

0 0

30

60

θ

90

[degree]

Fig. 3. Power of the emitted EM waves at the p- and 2p-cavity mode resonances in the [0 0 1], [1 0 1] and [0 0 1] directions.

N fying the Ohmic relation, I ¼ Nm bV with m = 0 (a), m = 1 (b) and m = 2 (c), where N is the number of junctions in the IJJ’s, i.e., N = 20 in the present case. It is seen that the I–V curve in this voltage region is split into three parts that nearly coincide with the Ohmic lines (a–c), respectively, that is, these three parts belong to different branches in the I–V characteristics on which the number of resistive junctions differs from one another. Note that the pmode resonance is situated in between the highest branch and the next highest one, as indicated by a red arrow in this figure. From this observation one understands that the strong emission stemming from the p-mode resonance takes place in the retrapping region where the highest voltage branch (uniform branch) becomes unstable, which is consistent with the experiments [3]. In conclusion our numerical simulations for the xz model well describe the THz wave emission due to the cavity-mode resonance in mesa samples of the Bi-2212 IJJ’s. The spatial pattern of the emitted EM waves and the I–V characteristics are consistent with

0.4

[100] direction

1.5

π−mode

Lx=0.5 λc α=0.1 η=104

π mode

internal Fiske step

S/Sp

I/jc

1

0.3

0.5

a

2 π mode

b c 0.2

0

0

5

10

15

20

V/Vp Fig. 2. Voltage dependence of the power of emitted EM waves in the [1 0 0] direction.

4

6

8

V/Vp Fig. 4. I–V characteristics in the low voltage region near the p-cavity mode resonance. The weak step structure seen in the high voltage branch is caused by the resonance with in-plane modes of spatial scale kJ  Lx (internal Fiske step).

T. Koyama et al. / Physica C 471 (2011) 1202–1205

experiments. The detail analysis of our numerical simulations and comparison with experiments will be published in a separate paper.

References [1] L. Ozyuzer, A.E. Koshelev, C. Kurter, N. Gopalsami, O.A. Li, M. Tachiki, K. Kadowaki, T. Yamashita, T. Tachiki, K.E. Gray, W.K. Kwok, U. Welp, Science 318 (2007) 291. [2] H. Minami, I. Kakeya, H. Yamaguchi, T. yamamoto, K. Kadowaki, Appl. Phys. Lett. 95 (2009) 232511.

1205

[3] M. Tsujimoto, K. Yamaki, K. Deguchi, T. Yamamoto, T. Kashiwagi, H. Minami, M. Tachiki, K. Kadowaki, Phys. Rev. Lett. 105 (2010) 037005. [4] H.B. Wang, S. Guenon, J. Yuan, A. Ishii, S. Arisawa, T. Hatano, T. Yamashita, D. Koelle, R. Kleiner, Phys. Rev. Lett. 102 (2009) 017006. [5] N. Bulaevskii, A.E. Koshelev, Phys. Rev. Lett. 99 (2007) 057002. [6] A.E. Koshelev, Phys. Rev. B 78 (2008) 174509. [7] S. Lin, X. Hu, Phys. Rev. Lett. 100 (2008) 247006. [8] X. Hu, S. Lin, Supercond. Sci. Technol. 23 (2010) 053001. [9] M. Tachiki, S. Fukuya, T. Koyama, Phys. Rev. Lett. 102 (2009) 127002. [10] T. Koyama, H. Matsumoto, M. Machida, K. Kadowaki, Phys. Rev. B 79 (2009) 104522. [11] Y. Nonomura, Phys. Rev. B 80 (2009) 140506. [12] N. Bulaevskii, A.E. Koshelev, J. Superecond. Novel Magn. 19 (2006) 349.