Transiently enhanced interlayer tunneling in optically driven high $T_c$ superconductors

TL;DR

This paper develops a theoretical model to explain the transient enhancement of interlayer tunneling in high-temperature superconductors induced by optical pulses, aligning with recent experimental observations and predicting larger effects under specific driving conditions.

Contribution

The paper introduces a pump-probe formalism for Josephson junction stacks to explain transient tunneling enhancement, highlighting the role of parametric driving near plasma frequencies.

Findings

Transient Josephson coupling can be enhanced with pulsed optical excitation.

Maximum effect occurs when driven above the plasma frequency.

Model reproduces experimental features below T_c.

Abstract

Recent pump-probe experiments reported an enhancement of superconducting transport along the $c$-axis of underdoped YBa$_2$Cu$_3$O$_{6+\delta}$ (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 is found experimentally. For bilayer Josephson junctions, an enhancement similar to that experimentally is predicted below the critical temperature $T_c$. This model reproduces the essential features of the enhancement measured below $T_c$. To reproduce the experimental results above $T_c$, we will explore extensions of this model, such as in-plane and amplitude fluctuations, elsewhere.

Equations

$$\ddot{\phi }+\gamma \dot{\phi }+{\omega }_{\text{Jp}}^2\left[1+A(t,t_e)\right]\sin \phi =I+\xi ,$$

$$A(t,t_e)=\frac{A_0}{2}\cos ({\omega }_{e}t+\varphi )\left[\tanh \left(\frac{t-t_e}{{\Delta }_e}\right)+1\right]$$

$$A(t,t_e)=A_0\cos ({\omega }_{e}t+\varphi )\exp \left[-\frac{(t-t_e{)}^2}{2{\Delta }_e^2}\right].$$

$$I(t-t_p)=I_0\cos \left[{\omega }_p(t-t_p)\right]\exp \left[-\frac{(t-t_p{)}^2}{2{\Delta }_p^2}\right],$$

$$V(t_s-t_p,t_s-t_e)={\int }_{-\infty }^{t_s}\rho (t_s-t^{\prime },t_s-t_e)I(t^{\prime }-t_p)d{t}^{\prime }.$$

$$V(\tau ,{\tau }_e)={\int }_{-\infty }^{\tau }\rho (\tau -t^{\prime },{\tau }_e)I(t^{\prime })d{t}^{\prime }.$$

$$\sigma (\omega ,{\tau }_e)\equiv \frac{1}{\rho (\omega ,{\tau }_e)}=\frac{I(\omega )d}{V(\omega ,{\tau }_e)}.$$

$$J_{\text{eff}}({\tau }_e)\equiv \frac{\hslash }{e^{\ast }d}\mathrm{Im}[\sigma (\omega ,{\tau }_e)\omega {]}_{\omega \to 0}.$$

$$J_{\text{eff}}^{\text{steady}}\simeq J_0\left[1+\frac{A_0^2{\omega }_{\text{Jp}}^2({\omega }_e^2-{\omega }_{\text{Jp}}^2)}{2({\omega }_e^2-{\omega }_{\text{Jp}}^2{)}^2+2{\gamma }^2{\omega }_e^2}\right].$$

$$\begin{array}{rl}A(t)& =A_0\cos ({\omega }_{e}t)+A_1\exp (i2{\omega }_{e}t)+A_1^{\ast }\exp (-i2{\omega }_{e}t)\\ & =A_0\cos ({\omega }_{e}t)+|A_1|\cos (2{\omega }_{e}t+{\varphi }_1),\end{array}$$

$$\begin{gathered} \left(\frac{-{\omega }^2-i\gamma \omega }{{\omega }_{\text{Jp}}^2}+1\right)\phi (\omega )+\frac{A_0}{2}\left[\phi (\omega +{\omega }_e)+\phi (\omega -{\omega }_e)\right] \\ +\frac{1}{2}\left[A_1\phi (\omega +2{\omega }_e)+A_1^{\ast }\phi (\omega -2{\omega }_e)\right]=\frac{I(\omega )}{{\omega }_{\text{Jp}}^2}. \end{gathered}$$

$$\frac{1}{2}\left[\begin{array}{cccccc}\ddots & & & & & \\ 2K_2& A_0& A_1^{\ast }& 0& 0& \\ A_0& 2K_1& A_0& A_1^{\ast }& 0& \\ A_1& A_0& 2K_0& A_0& A_1^{\ast }& \\ 0& A_1& A_0& 2K_{-1}& A_0& \\ 0& 0& A_1& A_0& 2K_{-2}& \\ & & & & \ddots & \end{array}\right]\left[\begin{array}{c}⋮\\ {\phi }_2\\ {\phi }_1\\ {\phi }_0\\ {\phi }_{-1}\\ {\phi }_{-2}\\ ⋮\end{array}\right]=\left[\begin{array}{c}⋮\\ 0\\ 0\\ \frac{I_0}{{\omega }_{\text{Jp}}^2}\\ 0\\ 0\\ ⋮\end{array}\right],$$

$$K_n\equiv \frac{-{\left({\omega }_p+n{\omega }_e\right)}^2-i\gamma ({\omega }_p+n{\omega }_e)}{{\omega }_{\text{Jp}}^2}+1.$$

$$\sigma ({\omega }_p)=i\left(\frac{e^{\ast }d}{\hslash }\right)\frac{I_0}{{\omega }_p\phi ({\omega }_p)}.$$

$$\frac{J_{\text{eff}}}{J_0}\simeq 1+A_0^2{\omega }_{\text{Jp}}^2\frac{2{\omega }_e^2+(|A_1|\cos {\varphi }_1-2){\omega }_{\text{Jp}}^2}{4({\omega }_e^2-{\omega }_{\text{Jp}}^2{)}^2-|A_1{|}^2{\omega }_{\text{Jp}}^4+4{\gamma }^2{\omega }_e^2},$$

$$\begin{gathered} \left[\begin{array}{c}\ddot{{\phi }_1}\\ \ddot{{\phi }_2}\end{array}\right]+\gamma \left[\begin{array}{c}\dot{{\phi }_1}\\ \dot{{\phi }_2}\end{array}\right]-\frac{4\pi e^{\ast }{\mu }^{2}I}{s}\left[\begin{array}{c}{\alpha }_1^{-1}\\ {\alpha }_2^{-1}\end{array}\right] \\ =\left[\begin{array}{cc}-(1+2{\alpha }_1){\Omega }_1^2& 2{\alpha }_2{\Omega }_2^2\\ 2{\alpha }_1{\Omega }_1^2& -(1+2{\alpha }_2){\Omega }_2^2\end{array}\right]\left[\begin{array}{c}{\phi }_1\\ {\phi }_2\end{array}\right]. \end{gathered}$$

$$\left(\frac{\hslash }{e^{\ast }}\right)\left[\begin{array}{c}\dot{{\phi }_1}\\ \dot{{\phi }_2}\end{array}\right]=\left[\begin{array}{cc}1+2{\alpha }_1& -2{\alpha }_2\\ -2{\alpha }_1& 1+2{\alpha }_2\end{array}\right]\left[\begin{array}{c}{V}_1\\ V_2\end{array}\right].$$

$$\sigma (\omega )=\frac{{ϵ}_{\text{av}}}{4\pi i}\frac{\left({\omega }^2+i\gamma \omega -{\omega }_{\text{Jp1}}^2\right)\left({\omega }^2+i\gamma \omega -{\omega }_{\text{Jp2}}^2\right)}{\omega \left({\omega }^2+i\gamma \omega -{\omega }_t^2\right)},$$

$$\bar{A}(t)=\frac{1}{2}\cos ({\omega }_{e}t+\varphi )\left[\tanh \left(\frac{t}{{\Delta }_e}\right)+1\right]e^{-{\gamma }_{e}t},$$

Full text

I Introduction

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.Först et al. (2015); Mankowsky et al. (2016); Nicoletti and Cavalleri (2016) Several examples are: optical switching of charge-density waves in transition metal dichalcogenides,Stojchevska et al. (2014) creation of effective magnetic fields in rare-earth compounds,Nova et al. (2016) and induction of lattice distortions in manganites.Rini et al. (2007); Tobey et al. (2008) In particular, in Refs. Fausti et al., 2011; Kaiser et al., 2014; Nicoletti et al., 2014; Först et al., 2014a, b; Hu et al., 2014; Mankowsky et al., 2014; Hunt et al., 2015; Casandruc et al., 2015; Mankowsky et al., 2015; Khanna et al., 2016; Hunt et al., 2016; Mankowsky et al., 2017; Hu et al., 2017, pump-probe techniques were used to control various layered high-$T_{c}$ superconductors. This resulted in the observations of light-enhanced and light-induced superconductivity. These intriguing experimental results were studied theoretically in Refs. Denny et al., 2015; Raines et al., 2015; Höppner et al., 2015; Patel and Eberlein, 2016; Okamoto et al., 2016; Sentef et al., 2017. 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 paper, we study the transient response of the superconducting phase below the critical temperature $T_{c}$ in layered systems, which we model as capacitively coupled Josephson junctions (see Fig. 1).van der Marel and Tsvetkov (1996); Koyama and Tachiki (1996); Matsumoto et al. (1999); Koyama (1999); Machida et al. (1999); Koyama (2000); van der Marel and Tsvetkov (2001); Koyama (2001, 2002); Machida and Koyama (2004); Machida and Sakai (2004); Shukrinov and Mahfouzi (2007); Shukrinov et al. (2009) The model is limited by its low dimensionality and lack of amplitude fluctuations of the order parameter, which prohibits us to describe light-induced superconductivity far above $T_{c}$. steady state properties of similar models have been investigated in Refs. Denny et al., 2015; Höppner et al., 2015; Okamoto et al., 2016. Here, in order to obtain the time-resolved conductivity, we introduce a pump-probe scheme similar to the one used experimentally by scanning through various pump-probe delay times with narrow probe pulses. In Sec. II, 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 also find that an effective critical temperature $T_{c}$ of the transient state, as defined below, can be larger than that of the steady state. In Sec. III, we first relate the transient behavior to driving the junction with additional higher harmonic frequencies. We then extend this analysis to point out an improved driving method that combines several harmonics in steady states. In Sec. IV, we use an effective model of a stack of weak and strong junctions, resembling the structure of YBCO.van der Marel and Tsvetkov (1996, 2001); Koyama (2002) Again, we find a transient enhancement of the Josephson coupling, and the comparison with experimental data shows qualitative agreement below $T_{c}$. Better quantitative description of the light-enhanced and -induced superconductivity needs to go beyond our model and include more complex physics such as amplitude fluctuations, lattice distortions, and competing charge order. Finally, Sec. V is the conclusion.

II Single Josephson junction: Transient dynamics

II.1 Model and Method

As our first model, we study a single Josephson junction with a bare Josephson coupling $J_{0}$, a thickness $d$, and a dielectric constant $\epsilon$. It has a characteristic plasma frequency $\omega_{\text{Jp}}=\sqrt{4\pi e^{*}dJ_{0}/\hbar\epsilon}$. The phase $\varphi$ of the junction obeys

[TABLE]

where $\gamma$ is a damping coefficient, $I$ an external current, and $\xi$ the thermal noise characterized by a temperature $T$ via $\langle\xi(t)\xi(t^{\prime})\rangle=2\gamma k_{B}T\delta(t-t^{\prime})$. We have included a parametric modulation of $J_{0}$ with an amplitude $A$ as $J_{0}\rightarrow J_{0}\left[1+A(t,t_{e})\right]$.Jung (1993); Zerbe et al. (1994); MacLachlan (1964); Nayfeh and Mook (2008); Citro et al. (2015); Zhu et al. (2016) As we will discuss in more detail later, modulation of $J_{0}$ may be induced by optically excited oxygen atoms inside the junction. Mathematically, the result does not change if the dielectric function $\epsilon$ or the interlayer thickness $d$ is modulated; they all periodically change $\omega_{\text{Jp}}$ and drive the junction parametrically. As the pump or driving pulse $A(t,t_{e})$, we choose either a continuous driving pulse with a nonzero rise time

[TABLE]

or a Gaussian pulse

[TABLE]

For both, $A_{0}$ is the amplitude of the driving, $\omega_{e}$ the driving frequency, $\phi$ the initial phase, and $\Delta_{e}$ the rise time or the pulse length, respectively. $t_{e}$ characterizes the starting time of the driving. The continuous driving gives access to the relaxation to the steady state, while the pulsed driving can illuminate short transient dynamics. We assume that the phase $\phi$ is uncontrolled, which is the case for the experiments discussed here.Jung (1993); Zerbe et al. (1994) In the following we always take a phase average over $\phi\in[0,2\pi]$.

In order to obtain a time-resolved conductivity, we follow the formulation of Ref. Kindt and Schmuttenmaer, 1999 (see also Refs. Němec et al., 2002; Orenstein and Dodge, 2015; Shao et al., 2016; Kennes et al., 2017a.) We add a probe pulse to the system,

[TABLE]

and then measure the voltage $V$ across the junction at sampling time $t_{s}$. We fix the pump time $t_{e}$ and scan $t_{p}$ and $t_{s}$. The number of probe pulses during a fixed time window and the shape of the spectrum determines the resolution of the obtained data. Without a driving pulse, the response of the system depends only on the difference $t_{s}-t_{p}$. However, with the time-dependent driving pulse, this is no longer the case (see Fig. 1), and the resistivity response $\rho$ becomes time dependent

[TABLE]

Moving to the relative time variables $\tau\equiv t_{s}-t_{p}$ and $\tau_{e}\equiv t_{s}-t_{e}$, we rewrite this as a convolution,

[TABLE]

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

[TABLE]

This quantity resembles the transient conductivity that was measured in Ref. Hu et al., 2014. As in Ref. Okamoto et al., 2016 we define an effective Josephson coupling $J_{\text{eff}}$ via

[TABLE]

This reduces to $J_{0}$ in equilibrium, and thus quantifies the effective interlayer tunneling energy.

II.2 Transient conductivity

In Fig. 2, we first show $\operatorname{Im}\sigma(\omega,\tau_{e})$ and $J_{\text{eff}}$ for continuous driving with $A_{0}=0.8$, $\Delta_{e}=15$, and $\gamma=0.1$ at $T=0$ (in the following, we put $\omega_{\text{Jp}}=1$). The probe pulse is taken as $I_{0}=0.1$, $\omega_{p}=0.1\omega_{\text{Jp}}$, and $\Delta_{p}=10$. We numerically integrate the equation of motion by Heun scheme with time step $h=10^{-3}$. As we have shown in Ref. Okamoto et al., 2016, the interlayer tunneling is enhanced (suppressed) at the blue- (red-) detuned side, and the driven steady state value is approximately

[TABLE]

Interestingly, for $\omega_{e}>\omega_{\text{Jp}}$, the transient value of $J_{\text{eff}}$ 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.

To elucidate the transient behavior further, we consider a Gaussian pulse, Eq. (3), with $A_{0}=0.8$ and $\omega_{e}=1.2\omega_{\text{Jp}}$. We plot $J_{\text{eff}}$ for several pump durations in Fig. 3(a) with $\gamma=0.1$. We find that the transient peak around $\tau_{e}=0$ is larger for smaller pump duration and that it is accompanied by a dip before and after the peak. When the duration is long enough, e.g., $\Delta_{e}=30$, we observe the enhancement of $J_{\text{eff}}$ only. Figure 3(b) shows the damping dependence of $J_{\text{eff}}$ at $\Delta_{e}=10$. When the damping is increased, the peak value decreases significantly, and at the same time the dip diminishes.

II.3 Temperature dependence

Figure 4(a) illustrates the temperature dependence of $J_{\text{eff}}$ for continuous driving starting at $\tau_{e}=0$. The parameters of the model are the same as the blue-detuned side of Fig. 2, and the results are averaged over $\sim 10^{5}$ samples. At low temperatures, after a few oscillations, the system approaches the driven steady state with a nonzero $J_{\text{eff}}$. As the temperature increases, while the peak still appears, it is shifted to earlier times. The following steady state has vanishingly small Josephson coupling (c.f. the purple line for $T=0.2$); the transient state has higher $T_{c}$ than the steady state. In Fig. 4(b), we compare $J_{\text{eff}}(T)$ of the undriven, steady state, and the transient (at $\tau_{e}=19$) values. We observe that the driven state has a lower $T_{c}$ compared to the undriven one. As we have shown in Ref. Okamoto et al., 2016, this is due to the fact that the fluctuations of the current $j(t)\equiv\sin\varphi(t)$ are increased when integrated over all frequencies $\int d\omega\langle|j(\omega)|^{2}\rangle$. This increased noise induces phase slips and destroys the phase coherence, even though the low-frequency part of the fluctuations, which determines $J_{\text{eff}}$, is reduced. 111Importance of phase slips was discussed recently in Ref. Kennes and Millis, 2017. A similar reduction of $T_{c}$ is seen in the bilayer case, when the plasma frequency of the weak junction is driven directly. Okamoto et al. (2016)

III Single Josephson junction: higher-order harmonic driving

In the previous section, we have demonstrated that the transient response of a single Josephson junction deviates from the steady state response under monochromatic driving. Here, we attribute the deviation to mixing of higher order driving frequencies, in particular $\pm n\omega_{e}$ with $n=2,3,\cdots$, that appear in the driving pulse. To elaborate on this observation, we discuss the response function of a Josephson junction that is driven by two frequencies. We present both a numerical calculation and an analytical estimate. These results expand on the response function that was derived in Ref. Okamoto et al., 2016 for a single driving frequency. We assume that the driving pulse has an admixture of the second harmonic frequency,

[TABLE]

with $\omega_{e}$ being the principal driving frequency. $A_{1}=|A_{1}|e^{i\phi_{1}}$ is the complex amplitude for the second order harmonic driving, where $\phi_{1}$ gives the phase difference between the first harmonic and the second one. By Fourier transforming and linearizing Eq. (1), we obtain

[TABLE]

We assume that the external current is monochromatic $I(\omega)=I_{0}\delta(\omega-\omega_{p})$ and the probing frequency $\omega_{p}$ is taken to be much smaller than the driving frequency $\omega_{e}$ and the plasma frequency $\omega_{\text{Jp}}$, since we are interested in the low frequency conductivity. This allows us to write down a discrete set of coupled equations for $\varphi_{n}\equiv\varphi(\omega_{p}+n\omega_{e})$ ($n\in\mathbb{Z}$) as

[TABLE]

where the diagonal elements are given by

[TABLE]

We obtain the solution by inverting the matrix numerically. We truncate the infinite matrix equation by taking $21$ modes ($n=-10,\dots,10$); we have checked that the convergence is well achieved in terms of the number of modes. Once we find $\varphi_{0}=\varphi(\omega_{p})$, we compute the conductivity from the Josephson relation $V=(\hbar/e^{*})\dot{\varphi}$ and $\sigma=Id/V$ as

[TABLE]

We define the effective Josephson coupling as in Eq. (8).

First, we discuss the case of in-phase driving, i.e., $\phi_{1}=0$. In Fig. 5, we plot $J_{\text{eff}}$ as a function of $|A_{1}|$ and $\omega_{e}$ for $A_{0}=0.2$, and $0.6$ at $\gamma=0.1$. We have excluded the regions where the driving pulse leads to Floquet parametric instability. Coddington and Carlson (1997); Cesari (2013) The stability is determined by the Floquet exponent obtained by integrating the equation of motion for one period of the driving, $2\pi/\omega_{e}$, with initial conditions $\left[\varphi(0),\dot{\varphi}(0)\right]=[1,0]$ and $[0,1]$.Coddington and Carlson (1997); Cesari (2013) We note that the instability region depends only weakly on the driving amplitude of the primal harmonic $A_{0}$. This indicates that the instability mainly comes from the second harmonic driving. Remarkably, the weak additional harmonic $A_{1}$ gives rise to larger values of $J_{\text{eff}}$ for the blue detuned side. As a competing effect, for larger values of $A_{1}$, the system reaches the primary Floquet instability lobe.

This result is modified by different values of the phase $\phi_{1}$. Figure 6 shows $J_{\text{eff}}$ as a function of $\omega_{e}$ at $A_{0}=0.4$, $|A_{1}|=0.1$, and $\gamma=0.1$ for various phase differences $\phi_{1}$. We find that the largest enhancement of $J_{\text{eff}}$ is obtained for in-phase driving $\phi_{1}=0$. As we dephase the two driving, the enhancement becomes weaker monotonically. To better understand this result, let us consider Eq.  (12) with only three modes ($n=-1,0,1$). The analytical expression of $J_{\text{eff}}$ is obtained as

[TABLE]

which is a generalization of Eq. (9). The second term is the correction by the driving. For $\omega_{e}>\omega_{\text{Jp}}$, the numerator becomes bigger when the two driving are near in-phase regime, $0\leq\phi_{1}<\pi/2$, while the near out-of-phase driving, $\pi/2\leq\phi_{1}\leq\pi$, leads to a smaller numerator. Overall tendency thus depends on the phase difference $\phi_{1}$. For small $A_{1}$, the numerator dominantly decides if $J_{\text{eff}}$ is enhanced or reduced compared to the $A_{1}=0$ case, since the denominator has only the quadratic contribution of $A_{1}$.

We note that these results suggest that $J_{\text{eff}}$ can be maximized in the driven steady state by carefully designing multifrequency 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., $\Delta_{e}=10$, the system is effectively driven by higher harmonics, in addition to the base frequency, which can lead to an initial suppression of $J_{\text{eff}}$, and then to strong increase of $J_{\text{eff}}$, as the higher harmonic admixture is reduced in time, crossing through the regime of optimal admixture.

IV Bilayer Josephson junctions

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 $d_{i}$, a dielectric constant $\epsilon_{i}$, a Josephson critical current $j_{i}$, and a bare plasma frequency $\Omega_{i}=\sqrt{4\pi e^{*}d_{i}j_{i}/\hbar\epsilon_{i}}$. We ignore fluctuations among different unit cells. Then the equation of motion of the phase differences $\varphi_{i}$ becomes,van der Marel and Tsvetkov (1996, 2001); Koyama (2002); Okamoto et al. (2016)

[TABLE]

$\alpha_{i}\equiv\epsilon_{i}\mu^{2}/sd_{i}$ is the capacitive coupling constant with $s$ being the thickness of the superconducting layer and $\mu$ the Thomas-Fermi screening length in the superconducting layers. The voltage is related to the phase differences by the generalized Josephson relations,Koyama and Tachiki (1996); Koyama (2002)

[TABLE]

For the undriven case at $T=0$, Eqs. (16) and (17) give

[TABLE]

where $\omega_{\text{Jp1, Jp2}}\simeq\Omega_{1,2}$ are the longitudinal plasma modes for weak and strong junctions, and $\omega_{t}\simeq\omega_{\text{Jp2}}$ is the transverse plasma mode.Koyama (2002) We take the parameters of the model as $\alpha_{1}=3$, $\alpha_{2}=1.5$, $\Omega_{1}=1$, $\Omega_{2}=12.5$, and $\gamma=0.1$. These are chosen to reproduce the ratio $\omega_{\text{Jp2}}/\omega_{\text{Jp1}}\sim 15.8$ of YBCO with appropriate $\alpha$ values for this compound of around $\sim 3$.Koyama (2002) We have $\omega_{\text{Jp1}}=1.58$ and $\omega_{\text{Jp2}}=25.1$. The probing pulse is taken as $I_{0}=0.1$, $\omega_{p}=1.5$, and $\Delta_{p}=1$ so that the frequencies around $\omega_{\text{Jp1}}$ are well resolved as the experimental condition.

In the experiment by Hu et al.,Hu et al. (2014) the optical pump pulse has a period of $\sim 50$fs and its duration is $\sim 0.3$ps, while the lifetime of the resonantly driven infrared $B_{1u}$ mode, which displaces apical oxygens along the $c$ axis, exceeds $2$ps.Subedi et al. (2014); Mankowsky et al. (2014, 2015) The oxygen motion primarily affects the interlayer motion of Cooper pairs, and thus we assume that the modulation of the Josephson critical currents $j_{i}$ derives from this driven phonon mode. The modulation of $j_{i}$ leads to parametric driving of the Josephson junction. As in the single Josephson junction case, we note that such parametric driving may be realized by other mechanisms such as modulation of the dielectric function,Denny et al. (2015). However, the microscopic origin of the modulation is not important for the following discussion. To imitate the transient phonon motions, we take the driving as

[TABLE]

with $\omega_{e}=26$, $\Delta_{e}=1$, and $\gamma_{e}=0.05$. This shows a sharp rise within several cycles of $\omega_{e}$ and then an exponential decay over dozens of cycles (Fig. 7). This parametric driving is included by changing the critical currents as $j_{1,2}\rightarrow j_{1,2}[1\pm a_{1,2}\bar{A}(t)]$; we assume that the driving is alternating along the junctions.

In Fig. 8, we compare the change of conductivity $\operatorname{Im}\Delta\sigma(\omega,\tau_{e})\equiv\operatorname{Im}\sigma(\omega,\tau_{e})-\operatorname{Im}\sigma^{\text{eq}}(\omega)$ obtained by simulations to the experimental result of Ref. Hu et al., 2014 at $T=10\text{K}\ll T_{c}$. For the simulation, we take $a_{1}=0.3$, $a_{2}=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 $\omega_{\text{Jp1}}$. We also observe that the transient change $\operatorname{Im}\Delta\sigma(\omega,\tau_{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 $\sim 0.5\omega_{\text{Jp1}}$ in the simulation is $\sim 15\%$ at most, while that in the experiment is $\sim 100\%$.222In equilibrium, $\operatorname{Im}\sigma(\omega=0.5\omega_{\text{Jp1}})\approx 10\ \Omega^{-1}\text{cm}^{-1}$. These discrepancies may arise due to physics that is not included in our simulation such as finite temperature effects, amplitude fluctuations of the order parameter, nonlinear lattice distortion,Mankowsky et al. (2014); Raines et al. (2015) and competing orders.Patel and Eberlein (2016); Sentef et al. (2017)

V Conclusions

In this paper, we have studied transient superconductivity in optically driven high $T_{c}$ superconductors using Josephson junction models below $T_{c}$. We find that 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 $T_{c}$, 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,Mankowsky et al. (2014); Raines et al. (2015) or competing charge order.Patel and Eberlein (2016); Sentef et al. (2017) 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. It is an interesting open question if the conductivity of phonon driven BCS superconductors Mitrano et al. (2016); Knap et al. (2016); Kim et al. (2016); Komnik and Thorwart (2016); Sentef et al. (2016); Kennes et al. (2017b); Murakami et al. (2017); Mazza and Georges (2017); Babadi et al. (2017); Nava et al. also shows larger enhancement when the higher-order harmonic driving is mixed.

Acknowledgements.

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.