Transiently enhanced interlayer tunneling in optically driven high ...

7 downloads 40 Views 6MB Size Report
Jun 14, 2017 - 4Department of Physics, Clarendon Laboratory, University of Oxford, Oxford OX1 3PU, United Kingdom. (Dated: June 15, 2017).
Transiently enhanced interlayer tunneling in optically driven high Tc superconductors Jun-ichi Okamoto,1, 2, ∗ Wanzheng Hu,3 Andrea Cavalleri,3, 4 and Ludwig Mathey1, 2

arXiv:1706.04554v1 [cond-mat.supr-con] 14 Jun 2017

1

Zentrum f¨ ur Optische Quantentechnologien and Institut f¨ ur Laserphysik, Universit¨ at Hamburg, 22761 Hamburg, Germany 2 The Hamburg Centre for Ultrafast Imaging, Luruper Chaussee 149, 22761 Hamburg, Germany 3 Max Planck Institute for the Structure and Dynamics of Matter, 22761 Hamburg, Germany 4 Department of Physics, Clarendon Laboratory, University of Oxford, Oxford OX1 3PU, United Kingdom (Dated: June 15, 2017) Recent pump-probe experiments reported an enhancement of superconducting transport along the c-axis of underdoped YBa2 Cu3 O6+δ (YBCO), induced by a mid-infrared optical pump pulse tuned to a specific lattice vibration. To understand this transient non-equilibrium state, we develop a pump-probe formalism for a stack of Josephson junctions, and we consider the tunneling strengths in presence of modulation with an ultrashort optical pulse. We demonstrate that a transient enhancement of the Josephson coupling can be obtained for pulsed excitation and that this can be even larger than in a continuously driven steady-state. Especially interesting is the conclusion that the effect is largest when the material is parametrically driven at a frequency immediately above the plasma frequency, in agreement with what found experimentally. For bilayer Josephson junctions, an enhancement similar to that experimentally is predicted below the critical temperature Tc . This model, which reproduces the essential features of the enhancement measured below Tc , does not, by design, reproduce the experimental results above Tc , a situation for which in-plane and amplitude superconducting fluctuations would be necessary to achieve a full quantitative understanding.

Recent pump-probe experiments have opened a new field in solid state physics by establishing a method to control material properties via laser pulses in the optical regime [1]. Several examples are: optical switching of charge-density waves in transition metal dichalcogenides [2], creation of effective magnetic fields in rareearth compounds [3], and induction of lattice distortions in manganites [4, 5]. In particular, in Refs. 6–19, pumpprobe techniques were used to control various layered high-Tc superconductors. This resulted in the observations of light-enhanced and light-induced superconductivity. These intriguing experimental results were studied theoretically in Refs. 20–25. However, these studies primarily focused on the steady-state of this driven system, while the experimental operation uses a pump pulse, with a pulse length that is typically around five times of the inverse optical frequency. It is therefore imperative to study the transient response of the driven system. In this work, we study the transient response of the superconducting phase below the critical temperature Tc in layered systems, which we model as coupled Josephson junctions (see Fig. 1). The model is limited by its low-dimensionality and lack of amplitude fluctuations of the order parameter, which prohibits us to describe lightinduced superconductivity far above Tc . In order to obtain the time-resolved conductivity, we use a pump-probe scheme similar to the one used experimentally by scanning through various pump-probe delay times with narrow probe pulses. We first consider a single Josephson junction as a simple model for the interlayer phase dynamics. When the frequency of the parametric driving is just above the Josephson plasma frequency, the effective Josephson coupling both in the transient and the driven steady-state is increased. In particular, when the driving pulse is narrow in time, the transient value can be larger

than the steady-state value. We explain this observation by considering the higher harmonics of the parametric driving. We also find that an effective critical temperature Tc of the transient state, as defined below, can be larger than that of the steady-state. Secondly, we use an effective model of a stack of weak and strong junctions, resembling the structure of YBCO [26, 27]. Again, we find a transient enhancement of the Josephson coupling, and the comparison with experimental data shows qualitative agreement below Tc . Better quantitative description of the light-enhanced and -induced superconductivity needs to go beyond our simple model and include more complex physics such as amplitude fluctuations, lattice distortions, and competing charge order. As our first model, we study a single Josephson junction with a bare Josephson coupling J0 , a thickness d, and a dielectric constant . It has a characteristic plasma frep quency ωJp = 4πe∗ dJ0 /~. The phase ϕ of the junction obeys 2 ϕ¨ + γ ϕ˙ + ωJp [1 + A(t, te )] sin ϕ = I + ξ,

(1)

where γ is a damping coefficient, I an external current, and ξ the thermal noise characterized by a temperature T via hξ(t)ξ(t0 )i = 2γkB T δ(t − t0 ). We have included a parametric modulation of J0 with an amplitude A as J0 → J0 [1 + A(t, te )] [28–33]. As the pump or driving pulse A(t, te ), we choose either a continuous driving pulse with a nonzero rise time     t − te A0 cos(ωe t + φ) tanh + 1 . (2) A(t, te ) = 2 ∆e or a Gaussian pulse   (t − te )2 A(t, te ) = A0 cos(ωe t + φ) exp − . 2∆2e

(3)

1 1 0

A(t)

2

A(t)

1 -10

A(t)-50 0 50 100 (b) 1 -10 0.1 A(t)-50 I(t-t )0 t p 100 = -50 50 p THz laser -10 0 CuO2 0.1 t p 100 = -50 50 V(t) -50 I(t-t p )0 -0.1 -10 0.1 t = -50 -50 0 50 100 I(t-t ) 50 V(t) -50 p 100 p0 -0.1 0 0.1 =0 -50 CuO2 0.1 V(t) -50 I(t-t )0 tt pp 100 = 50 p -0.1 0 0 0.1 t p 100 =0 0 50 V(t) -50 -0.1 -0.1 0 0.1 t =0 -50 0 50 100 -50 0 50 p 100 -0.1 0 0.1 t = 0.1 t p 100 =0 50 -50 0 50 p -0.1 apical oxygen 0.100 t 100 = 50 -50 0 50 p -0.1 -0.1 0 0.1 t 100 = 50 -50 0 50 100 -50 0 50 p -0.1 0 0.1 t t = 50 -50 0 50 p 100 -0.1 0 t -50 0 50 100 FIG. 1. (a) Schematic -0.1 depiction of YBCO. Superconducting t -50 of bilayer 0 50 100 CuO2 layers (grey) form a stack Josephson junc-

(a)

We add a probe pulse to the system,   (t − tp )2 , I(t − tp ) = I0 cos [ωp (t − tp )] exp − 2∆2p

and then measure the voltage V across the junction at sampling time ts . We fix the pump time te and scan tp and ts . Without a driving pulse, the response of the system depends only on the difference ts − tp . However, with the time-dependent driving pulse, this is no longer the case (See Fig. 1), and the resistivity response ρ becomes time-dependent Z ts V (ts −tp , ts −te ) = ρ(ts −t0 , ts −te )I(t0 −tp )dt0 . (5) −∞

t tions. THz pulses (wavy lines) excite apical oxygen atoms (circles) that induce oscillations of j1 and j2 . (b) Typical time-dependent voltage response V (t) (solid lines) for different probe pulses I(t − tp ) (dashed lines) at tp = −50, 0 and 50. The driving amplitude A(t) is also depicted.

Moving to the relative time variables τ ≡ ts − tp , and τe ≡ ts − te , we rewrite this as a convolution, Z τ V (τ, τe ) = ρ(τ − t0 , τe )I(t0 )dt0 . (6) −∞

10 5

Fourier transforming the above equation in terms of τ , we define the time-dependent conductivity as

0

0.8 1

σ(ω, τe ) ≡

-1

/ Jp1

(b1) Imσ, ωe = 0.8ωJp

cm -1

0.6 (a1) Imσ, ωe = 1.2ωJp

(4)

I(ω) 1 = . ρ(ω, τe ) V (ω, τe )

(7)

-5

This quantity resembles the transient conductivity that -10 was measured in Ref. 11. As in Ref. 24 we define an effective Josephson coupling Jeff via 1

0

1

2

3

J eff / J 0

0.8

1.2 (b) theory

0.5 0

2

80

e

3

0 -20

/2 e Jp1

4

0

20

5

40

60

6

80

e

0

0

0.5

20

20

FIG. 2. Transient imaginary conductivity Im σ(ω, τe ) and Jeff5(τ(a) driving at 5T(b) = 0. (a) Blue-detuned 0 e ) for continuous = 10 = 0.1 e case4 ωe = 1.2ωJp . (b) Red-detuned 4 case ωe = 0.8ωJp . = 0.2 = 20 J eff / J 0

2 1

60

(b) -1

1

= 0.5

3

40

= 30

60

2

e

10 0 -10

-0.5

40

J eff / J 0

e

e

3

80

80

0 20 of the 40 driving, 60 80 For0 both, A0 is1-20 the ωe the drivA(t)amplitude 0 0 -40 -20 20 40 -20 0 20 time or 40 ing frequency, φ the initial phase, and ∆e the rise e the startthe pulse length, -1e respectively. te characterizes -100 1 -20 0 20 40 60 80 5 10 ing time (c)J of effthe The continuous driving gives (A1,driving. ωe) e) V(τ,τ 0.25 access to the-100 relaxation to the steady-state, while40.51the -50 5 10 unstable pulsed0.2driving can illuminate short transient dynamics. 0 0 3 0 0.15 -100 We assume that0 the phase φ is uncontrolled, which is 0the 100 -1 2 -100 case for discussed the -20 0 20 40 here 60 [28,8029]. In 0.1 the experiments -5 50 -0.5 following we always take a phase 2π]. e average over φ ∈ [0, 1 0.05 In order to obtain a time-resolved conductivity, we fol-10 0 100 low the formulation of0Ref. 2034 also Refs. 35 and0 -1 36.) 0.8 -20 0.9 1 (see 1.280 40 1.1 60 -100

1

0

I(t)

-1 -100

A1

0

τ

-50 50 100

e

J0

2 1.5

T=0 T=0.1 T=0.15

τe Jpe

2.5

(a)

2 0

2.5

/

(8)

(b)

Undriven Steady state Transient peak

1

60

~ Im[σ(ω, τe )ω]ω=0 . e∗ d

This reduces to J0 in the equilibrium, and thus quantifies effective interlayer tunneling energy. In Fig. 2, we first show Im σ(ω, τe ) and Jeff for continuous driving with A0 = 0.8, ∆e = 15, and γ = 0.1 at T = 0 (In the following, we put ωJp = 1). The probe pulse is taken as I0 = 0.1, ωp = 0.1ωJp , and ∆p = 10. We numerically integrate the equation of motion by Heun scheme with time step h = 10−3 . As we have shown in Ref. 24, the interlayer tunneling is enhanced (supressed) at the blue- (red-) detuned side, and the driven steadystate value is approximately " # 2 2 A20 ωJp (ωJp − ωe2 ) steady Jeff ' J0 1 − (9) 2 − ω 2 )2 + 2γ 2 ω 2 . 2(ωJp e e 0

140

-1 -20

1

20

-1 0

0

1

-100

0

Jeff (τe ) ≡

-0.5 the

100 -20

0 -20

6

2

1 1

5

(b2) Jeff, ωe = 0.8ωJp

J eff / J 0

/ Jp1

(a2)0.6 Jeff, ωe =1.2 ωJp 2

4

/2 e Jp1

cm -1

0.4

-1

1.2 (a) exp.

Interestingly, for ωe > ωJp , the transient value of Jeff first shows a dip in the initial stage of the driving, followed by a large peak, and then reaches to the steady-state after a few small oscillations. We will explain this behavior in more detail below. -50 To elucidate 0 50the transient 100 behavior further, we consider a Gaussian pulse, Eq. (3), with A0 = 0.8 and ωe = 1.2ωJp . We plot Jeff for several pump durations in Fig. 03(a) with γ = 0.1. We find that the transient -50 50 100

20

40

60

0 -20

80

0

20

40

e

(a)

J eff / J 0

4

0

80

0.8

e e

3

e

5

= 10

4

= 20 = 30

2 1 -20

(b)

= 0.1 = 0.2 = 0.5

3 2

0

20

0 -20

40

0

20

40

e

A1

1.5

0

1.2

Jp

3 2.5

(a)

T=0 T=0.1 T=0.15 T=0.2

(b)

Undriven Steady state Transient

2

1

0 -10

/

1.1

1.5 1 0.5

0

10

20

30

40

0 0

0.1

0.2

0.3

0.4

Temperature

5 4

unstable 3

0.15

2

0.1

FIG. 4. (a) Temperature dependence of Jeff (τe ) for continuous driving. We use γ = 0.1, ∆e = 15, and ωe = 1.2. (b) Temperature dependence of Jeff for A0 = 0 (undriven), A0 = 0.8 (steady-state), and A0 = 0.8 (transient value at τe = 19).

1

0.05 0 0.8

0.9

1 e

/

1.1

0

1.2

Jp

J eff / J 0

J eff / J 0

FIG. 3. (a) Jeff (τe ) for several pump widths ∆e for γ = 0.1 2.5 damping factors γ for at 2.5 T = T=0 0. (b) Jeff (τe ) for (a) different (b) Undriven Steady state ∆e =2 10 T=0.1 at T = 0. (c) Jeff /J0 in the2 driven steady-state with T=0.15 Transient peak additional second harmonic driving A cos(2ω t) as a function 1.5 1 e 1.5 T=0.2 of A11 and ωe for A0 = 0.6 and γ = 0.1. 0.5 0 -10

2

e

(c)Jeff (A1, ωe)

0.2

2.5

0.5

e

0.25

1 e

1

0 -40

0.9

e

J eff / J 0

5

60

J eff / J 0

0

J eff / J 0

0 -20

1

0.5 0

10

20

30

40

0 0

0.1

0.2

0.3

0.4

e 0 is larger for smaller Temperature peak around τe = pump duration, and that it is accompanied by a dip before and after the peak. When the duration is long enough, e.g. ∆e = 30, we observe the enhancement of Jeff only. Fig. 3(b) shows the damping dependence of Jeff at ∆e = 10. When the damping is increased, the peak value decreases significantly, and at the same time the dip diminishes. We attribute this transient behavior to mixing in higher order driving frequencies, in particular ±nωe with n = 2, 3, · · · . As an example, we have performed a steady-state calculation within the linear approximation that was used in Ref. 24, with an additional second harmonic driving term A1 cos(2ωe t) (for the details, see Supplemental Material.) The resulting Jeff is shown in Fig. 3(c), as a function of A1 and ωe , and for fixed A0 = 0.6 and γ = 0.1. We have excluded the dynamically unstable regions determined by Floquet analysis [37, 38]. Remarkably, the weak additional harmonic A1 gives rise to larger values of Jeff for the blue detuned side. As a competing effect, for larger values of A1 , the system reaches the primary Floquet instability lobe [37, 38]. We note that this result suggests that Jeff can be maximized in the driven steady-state by carefully designing multi-frequency optical driving, which takes advantage of the increase that can be achieved by adding higher harmonics, while avoiding the parametric instability regime. These two features also compete in the transient response due to a short driving pulse. In particular, at the initial stage of a short pulse, e.g. ∆e = 10, the system is effectively driven by higher harmonics, in addition to the base frequency, which can lead to an initial suppression

of Jeff , and then to strong increase of Jeff , as the higher harmonic admixture is reduced in time, crossing through the regime of optimal admixture. Fig. 4(a) illustrates the temperature dependence of Jeff for continuous driving. The parameters are the same as the blue-detuned side of Fig. 2, and the results are averaged over ∼ 105 samples. At low temperatures, after a few oscillations, the system approaches the driven steady-state with a nonzero Jeff . As the temperature increases, while the peak still appears (shifted to an earlier time), the corresponding steady-state has vanishingly small Josephson coupling (c.f. the purple line for T = 0.2); the transient state has higher Tc than the steady-state. In Fig. 4(b), we compare Jeff (T ) of the undriven, steady-state, and the transient (at τe = 19) values. We observe that the driven state has a lower Tc compared to the undriven one. As we have shown in Ref. 24, this is due to the fact that the fluctuations of the current j(t) ≡ sinRϕ(t) are increased when integrated over all frequencies dωh|j(ω)|2 i. This induces phase slips and destroys the phase coherence, even though the low-frequency part of the fluctuations, which determines Jeff , is reduced. A similar reduction of Tc is seen in the bilayer case, when the plasma frequency of the weak junction is directly driven [24]. As our second model, we consider a stack of alternating weak (i = 1) and strong (i = 2) junctions, see Fig. 1. Each junction is characterized by a thickness di , a dielectric constant i , a Josephsonp critical current ji , and a bare plasma frequency Ωi = 4πe∗ di ji /~i . We ignore fluctuations among different unit cells. Then the equation of motion of the phase differences ϕi becomes [24, 26, 27],       4πe∗ µ2 I α1−1 ϕ¨1 ϕ˙1 +γ − ϕ¨2 ϕ˙2 α2−1 s    −(1 + 2α1 )Ω21 2α2 Ω22 ϕ1 = . (10) 2α1 Ω21 −(1 + 2α2 )Ω22 ϕ2 αi ≡ i µ2 /sdi is the capacitive coupling constant with s being the thickness of the superconducting layer and µ

4

5

0.9

0

1 1.1 (a) exp.

-5

0.8

0

1

-0.5

1.2 (b) theory 0

cm -1

0.5

2

4 Jp1

-1

8

/2 1

0

-1 -20

0

1

-100

e

6

0

FIG. 5. Changes in transient imaginary conductivity Im ∆σ(ω, τe ).0.5 The simulated conductivity is rescaled to fit the value at ω = 0.5ωJp1 in equilibrium to the experimental 0 one. (a) Experimental data from Ref. 11 at T = 10 K  Tc . (b) A simulated result at γ = 0.1, a1 = 0.3, and a2 = 0.6 at -0.5 T = 0. 20

e

20

40

60

80

0

80

A(t)

60

60

1-20

10 0 -10

(b) -1

40

40

80

0

1

0

-1 -100

I(t)

-1 1 that in the -100 experiment is [41]. These discrepan-20 0 20 ∼ 100% 40 60 80 10 V(τ,τe) cies may arise due to physics that is not included in our -50 0.5 -100 1 5 10 simulation such as finite temperature effects, amplitude 0 0 0 fluctuations of0 the order parameter, dis-100 nonlinear lattice 0 100 -1 tortion [12, 21], [23, 25].-100 -20 and 0 competing 20 40 orders 60 80 -50 -5

50

50

0

τ

-50

0

50

-0.5

e In conclusion, we have studied transient superconduc-10 tivity in optically driven high T usc superconductors 100 -1 0 20 40 60 80 -100 -50 ing Josephson-20junction models below T . We find that c τe e the transient state shows enhanced interlayer tunneling, which can be larger than the steady-state value, when the system is driven near the blue-detuned side of the Josephson plasma frequency. We have explained the transient behavior by considering the higher order harmonics in driving. We have also shown that our bilayer model can phenomenologically explain the temporal change of the imaginary part of conductivity seen in experiments on YBCO below Tc , while quantitative differences still remain; in particular it can hardly explain the light-induced superconductivity. The differences may derive from more complex physics including amplitude fluctuations, lattice distortion [12, 21], or competing charge order [23, 25]. We have also demonstrated that admixing higher harmonics in the driving operation can result in an additional enhancement of the c-axis transport. This observation opens the door towards optimal control of superconductivity via optical driving, by combining several higher harmonics.

100

with ωe = 26, ∆e = 1, and γe = 0.05. This shows a sharp rise within several cycles of ωe , and then an exponential decay over dozens of cycles. This parametric driving is included by changing the critical currents as ¯ j1,2 → j1,2 [1 ± a1,2 A(t)]; we assume that the driving is alternating along the junctions. In Fig. 5, we compare the change of conductivity Im ∆σ(ω, τe ) ≡ Im σ(ω, τe ) − Im σ eq (ω) obtained by simulations to the experimental result of Ref. 11 at T = 10K  Tc . For the simulation, we take a1 = 0.3, a2 = 0.6, and T = 0. At low frequencies, on the rise of the driving, a peak appears after a small dip in the simulation. This is similar to the single junction case. The peak is followed by a decay over few oscillations, relaxing back to the original state. The period of such oscillations is approximately one cycle of ωJp1 . We also observe that the transient change Im ∆σ(ω, τe ) becomes negative at high frequencies in the simulation. This overall transient behavior of the simulation is qualitatively similar to the experimental one, while a few discrepancies remain. For instance, a dip in the initial stage of the driving and subsequent small oscillations are absent in the experiment. Also, the relative enhancement at low frequencies ∼ 0.5ωJp1 in the simulation is ∼ 15% at most, while

/2

20

(13)

Jp1

-10 1

8

0

    t 1 ¯ + 1 e−γe t , A(t) = cos(ωe t + φ) tanh 2 ∆e

e

0.6

6

100 -20

where ωJp1, Jp2 ' Ω1,2 are the longitudinal plasma modes for weak and strong junctions, and ωt ' ωJp2 is the transverse plasma mode [27]. We take the parameters of the model as α1 = 3, α2 = 1.5, Ω1 = 1, Ω2 = 12.5, and γ = 0.1. These are chosen to reproduce the ratio ωJp2 /ωJp1 ∼ 15.8 of YBCO with appropriate α values for this compound of around ∼ 3 [27]. We have ωJp1 = 1.58 and ωJp2 = 25.1. The probing pulse is taken as I0 = 0.1, ωp = 1.5, and ∆p = 1 so that the frequencies around ωJp1 are well resolved as the experimental condition. In the experiment by Hu et al. [11], the optical pump pulse has a period of ∼ 50fs and its duration is ∼ 0.3ps, while the lifetime of the resonantly driven infrared B1u mode, which displaces apical oxygens along the c-axis, exceeds 2ps [12, 15, 40]. We assume that the parametric driving derives from this driven phonon mode. To imitate the transient phonon motions, we take the driving as

/ Jp1

For the undriven case at T = 0, Eqs. (10) and (11) give  2  2 2 ω + iγω − ωJp2 av ω 2 + iγω − ωJp1 , (12) σ(ω) = 4πi ω (ω 2 + iγω − ωt2 )

4

-1

2

0.4

-1

0.8

cm -1

10

0.7

/ Jp1

the Thomas-Fermi screening length in the superconducting layers. The voltage is related to the phase differences by the generalized Josephson relations [27, 39],       ~ ϕ˙1 1 + 2α1 −2α2 V1 = . (11) ∗ ϕ ˙ −2α 1 + 2α V2 e 2 1 2

We acknowledge the support from the Deutsche Forschungsgemeinschaft (through SFB 925 and EXC 1074) and from the Landesexzellenzinitiative Hamburg, which is supported by the Joachim Herz Stiftung. J.O. thanks S. A. Sato for helpful discussions about details of the pump-probe formalism.

0

50

5



[email protected] [1] R. Mankowsky, M. F¨ orst, and A. Cavalleri, Reports Prog. Phys. 79, 064503 (2016). [2] L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S. Brazovskii, and D. Mihailovic, Science 344, 177 (2014). [3] T. F. Nova, A. Cartella, A. Cantaluppi, M. F¨ orst, D. Bossini, R. V. Mikhaylovskiy, A. V. Kimel, R. Merlin, and A. Cavalleri, Nat. Phys. 13, 132 (2016). [4] M. Rini, R. Tobey, N. Dean, J. Itatani, Y. Tomioka, Y. Tokura, R. W. Schoenlein, and A. Cavalleri, Nature 449, 72 (2007). [5] R. I. Tobey, D. Prabhakaran, A. T. Boothroyd, and A. Cavalleri, Phys. Rev. Lett. 101, 197404 (2008). [6] D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Science 331, 189 (2011). [7] S. Kaiser, C. R. Hunt, D. Nicoletti, W. Hu, I. Gierz, H. Y. Liu, M. Le Tacon, T. Loew, D. Haug, B. Keimer, and A. Cavalleri, Phys. Rev. B 89, 184516 (2014). [8] D. Nicoletti, E. Casandruc, Y. Laplace, V. Khanna, C. R. Hunt, S. Kaiser, S. S. Dhesi, G. D. Gu, J. P. Hill, and A. Cavalleri, Phys. Rev. B 90, 100503 (2014). [9] M. F¨ orst, A. Frano, S. Kaiser, R. Mankowsky, C. R. Hunt, J. J. Turner, G. L. Dakovski, M. P. Minitti, J. Robinson, T. Loew, M. Le Tacon, B. Keimer, J. P. Hill, A. Cavalleri, and S. S. Dhesi, Phys. Rev. B 90, 184514 (2014). [10] M. F¨ orst, R. I. Tobey, H. Bromberger, S. B. Wilkins, V. Khanna, A. D. Caviglia, Y. D. Chuang, W. S. Lee, W. F. Schlotter, J. J. Turner, M. P. Minitti, O. Krupin, Z. J. Xu, J. S. Wen, G. D. Gu, S. S. Dhesi, A. Cavalleri, and J. P. Hill, Phys. Rev. Lett. 112, 157002 (2014). [11] W. Hu, S. Kaiser, D. Nicoletti, C. R. Hunt, I. Gierz, M. C. Hoffmann, M. Le Tacon, T. Loew, B. Keimer, and A. Cavalleri, Nat. Mater. 13, 705 (2014). [12] R. Mankowsky, A. Subedi, M. F¨ orst, S. O. Mariager, M. Chollet, H. T. Lemke, J. S. Robinson, J. M. Glownia, M. P. Minitti, A. Frano, M. Fechner, N. A. Spaldin, T. Loew, B. Keimer, A. Georges, and A. Cavalleri, Nature 516, 71 (2014). [13] C. R. Hunt, D. Nicoletti, S. Kaiser, T. Takayama, H. Takagi, and A. Cavalleri, Phys. Rev. B 91, 020505 (2015). [14] E. Casandruc, D. Nicoletti, S. Rajasekaran, Y. Laplace, V. Khanna, G. D. Gu, J. P. Hill, and A. Cavalleri, Phys. Rev. B 91, 174502 (2015). [15] R. Mankowsky, M. F¨ orst, T. Loew, J. Porras, B. Keimer, and A. Cavalleri, Phys. Rev. B 91, 094308 (2015). [16] V. Khanna, R. Mankowsky, M. Petrich, H. Bromberger, S. A. Cavill, E. M¨ ohr-Vorobeva, D. Nicoletti, Y. Laplace, G. D. Gu, J. P. Hill, M. F¨ orst, A. Cavalleri, and S. S. Dhesi, Phys. Rev. B 93, 224522 (2016). [17] C. R. Hunt, D. Nikoletti, S. Kaiser, D. Pr¨ opper, T. Loew, J. Porras, B. Keimer, and A. Cavalleri, Phys. Rev. B 94, 224303 (2016). [18] R. Mankowsky, M. Fechner, M. F¨ orst, A. von Hoegen, J. Porras, T. Loew, G. L. Dakovski, M. Seaberg, S. M¨ oller, G. Coslovich, B. Keimer, S. S. Dhesi, and A. Cavalleri, Struct. Dyn. 4, 044007 (2017). [19] W. Hu, D. Nicoletti, A. V. Boris, B. Keimer, and A. Cavalleri, Phys. Rev. B 95, 104508 (2017).

[20] S. J. Denny, S. R. Clark, Y. Laplace, A. Cavalleri, and D. Jaksch, Phys. Rev. Lett. 114, 137001 (2015). [21] Z. M. Raines, V. Stanev, and V. M. Galitski, Phys. Rev. B 91, 184506 (2015). [22] R. H¨ oppner, B. Zhu, T. Rexin, A. Cavalleri, and L. Mathey, Phys. Rev. B 91, 104507 (2015). [23] A. A. Patel and A. Eberlein, Phys. Rev. B 93, 195139 (2016). [24] J. I. Okamoto, A. Cavalleri, and L. Mathey, Phys. Rev. Lett. 117, 227001 (2016). [25] M. A. Sentef, A. Tokuno, A. Georges, and C. Kollath, Phys. Rev. Lett. 118, 087002 (2017). [26] D. van der Marel and A. A. Tsvetkov, Phys. Rev. B 64, 024530 (2001). [27] T. Koyama, J. Phys. Soc. Japan 71, 2986 (2002). [28] P. Jung, Phys. Rep. 234, 175 (1993). [29] C. Zerbe, P. Jung, and P. H¨ anggi, Phys. Rev. E 49, 3626 (1994). [30] N. W. MacLachlan, Theory and Application of Mathieu Functions (Dover, 1964). [31] A. H. Nayfeh and D. T. Mook, Nonlinear Oscillations (Wiley, 2008). [32] R. Citro, E. G. Dalla Torre, L. D’Alessio, A. Polkovnikov, M. Babadi, T. Oka, and E. Demler, Ann. Phys. (N. Y). 360, 694 (2015). [33] B. Zhu, T. Rexin, and L. Mathey, Zeitschrift fur Naturforsch. - Sect. A J. Phys. Sci. 71, 921 (2016). [34] J. T. Kindt and C. A. Schmuttenmaer, J. Chem. Phys. 110, 8589 (1999). [35] H. Nˇemec, F. Kadlec, and P. Kuˇzel, J. Chem. Phys. 117, 8454 (2002). [36] J. Orenstein and J. S. Dodge, Phys. Rev. B 92, 134507 (2015). [37] A. Coddington and R. Carlson, Linear Ordinary Differential Equations (Society for Industrial and Applied Mathematics, 1997). [38] L. Cesari, Asymptotic Behavior and Stability Problems in Ordinary Differential Equations (Springer Berlin Heidelberg, 2013). [39] T. Koyama and M. Tachiki, Phys. Rev. B 54, 16183 (1996). [40] A. Subedi, A. Cavalleri, and A. Georges, Phys. Rev. B 89, 220301 (2014). [41] In equilibrium, Im σ(ω = 0.5ωJp1 ) ≈ 10 Ω−1 cm−1 .

1

Supplemental Material: Transiently enhanced interlayer tunneling in optically driven high Tc superconductors STEADY STATE CONDUCTIVITY UNDER DRIVING WITH HIGHER HARMONICS

Here we present the details of the steady state calculations of conductivity under periodic driving with higher harmonics. We start from a sine-Gordon equation with a parametric driving pulse, 2 ϕ¨ + γ ϕ˙ + ωJp [1 + A(t)] sin ϕ = I,

(S1)

where ϕ is the superconducting phase difference of the single Josephson junction, ωJp the plasma frequency of the junction, γ a damping coefficient, and I an external current. As a driving pulse, we consider the admixture of two harmonics, A(t) = A0 cos(ωe t) + A1 cos(2ωe t + φ1 ),

(S2)

with ωe being the principal driving frequency. The second term is the second order harmonic, and φ1 is the phase difference between the first harmonic and the second one. We find that the results do not depend on the phase difference φ1 , so we will omit φ1 in the following. By Fourier transforming and linearizing Eq. (S1), we obtain ! A0 A1 I(ω) −ω 2 − iγω + 1 ϕ(ω) + [ϕ(ω + ωe ) + ϕ(ω − ωe )] + [ϕ(ω + 2ωe ) + ϕ(ω − 2ωe )] = 2 . (S3) 2 ωJp 2 2 ωJp We assume that the external current is monochromatic I(ω) = I0 δ(ω − ωp ) and the probing frequency ωp is taken to be much smaller than the driving frequency ωe and the plasma frequency ωJp , since we are interested in the low frequency conductivity. This allows us to write down a discrete set of coupled equations for ϕn ≡ ϕ(ωp + nωe ) (n ∈ Z) as      .. .. .. . . .          A0 2K2 A0 A1 0 0 0   ϕ2   0     ϕ1   0  A1 A0 2K1 A0 A 0 0 1     1 2   ϕ0  = I0 /ωJp  0 A1 A0 2K0 A0 A1 0 (S4)  ,     2  ϕ−1   0  0 0 A A 2K A A 1 0 −1 0 1          0 0 0 A1 A0 2K−2 A0   ϕ−2   0   .. .. .. . . . where the diagonal elements are given by Kn ≡

! −ωp2 − iγωp +1 . 2 ωJp

(S5)

5

0.25

5

A0 = 0.2

0.25

A0 = 0.6

4

4

0.2 3

0.15

A1

A1

0.2

0.1

2

unstable

1

0.05

0

0

0

0.9

1 e

/

1.1 Jp

1.2

unstable

0.1

0.05 0.8

3

0.15

2 1

0.8

0.9

1 e

/

1.1

1.2

0

Jp

FIG. S1. Jeff /J0 in the driven steady-state with additional second harmonic driving A1 cos(2ωe t) as functions of A1 and ωe for γ = 0.1 and A0 = 0.2 (left) or 0.6 (right).

2 The solution can be easily obtained by numerically inverting the matrix. In practice, we truncate the infinite matrix equation by taking 21 modes (n = −10, . . . , 10); we have checked that the convergence is well achieved in terms of the number of modes. Once we find ϕ0 = ϕ(ωp ), we can compute the conductivity from the Josephson relation V = (~/e∗ )ϕ˙ and σ = Id/V (d: the thickness of the junction) as  ∗  I0 e d . (S6) σ(ωp ) = i ~ ωp ϕ(ωp ) We define the effective Josephson coupling as Jeff

  ~ iI0 ≡ ∗ Im[σ(ω)ω]ω=0 = Im . e d ϕ(ωp → 0)

(S7)

In Fig. S1, we plot the obtained Jeff as functions of A1 and ωe for A0 = 0.2, and 0.6 at γ = 0.1. We have excluded the regions where the driving pulse leads to Floquet parametric instability. The stability is determined by the Floquet exponent obtained by integrating the equation of motion for one period of the driving, 2π/ωe , with initial conditions [ϕ(0), ϕ(0)] ˙ = [1, 0] and [0, 1] [37, 38]. We note that the instability region does not depend much on the driving amplitude of the primal harmonic A0 . This indicates that the instability mainly comes from the second harmonic driving. In principle, we can further generalize the analysis including even higher harmonics. This may result in larger enhancement of Jeff , while, as we have seen, it is important to avoid the unnecessary heating caused by the Floquet instability.