Evidence for time-reversal symmetry breaking of the superconducting state near twin-boundary interfaces in FeSe

Junctions and interfaces consisting of unconventional superconductors provide an excellent experimental playground to study exotic phenomena related to the phase of the order parameter. Not only the complex structure of unconventional order parameters have an impact on the Josephson effects, but also may profoundly alter the quasi-particle excitation spectrum near a junction. Here, by using spectroscopic-imaging scanning tunneling microscopy, we visualize the spatial evolution of the local density of states (LDOS) near twin boundaries (TBs) of the nodal superconductor FeSe. The $\pi/2$ rotation of the crystallographic orientation across the TB twists the structure of the unconventional order parameter, which may, in principle, bring about a zero-energy LDOS peak at the TB. The LDOS at the TB observed in our study, in contrast, does not exhibit any signature of a zero-energy peak and an apparent gap amplitude remains finite all the way across the TB. The low-energy quasiparticle excitations associated with the gap nodes are affected by the TB over a distance more than an order of magnitude larger than the coherence length $\xi_{ab}$. The modification of the low-energy states is even more prominent in the region between two neighboring TBs separated by a distance $\approx7\xi_{ab}$. In this region the spectral weight near the Fermi level ($\approx\pm$0.2~meV) due to the nodal quasiparticle spectrum is almost completely removed. These behaviors suggest that the TB induces a fully-gapped state, invoking a possible twist of the order parameter structure which breaks time-reversal symmetry.


I. INTRODUCTION
When two superconductors are in close proximity, they are influenced by each other via the tunneling of Cooper pairs. The Cooper-pair tunneling results in the flow of a superconducting Josephson current, which has been studied for decades and is used in various superconducting quantum devices [1]. The Josephson current is governed by the phase difference of the order parameters of the two superconductors. Therefore, Josephson junctions consisting of unconventional superconductors, where the superconducting order parameter changes its sign depending on the momentum direction, serve as a unique platform where novel phase-related phenomena, e.g., spontaneous formation of half flux quanta in a tri-junction of cuprate superconductors [2], take place. Compared with the well-investigated Josephson currents, the spatial and energy dependence of the superconducting order parameter and quasiparticle states around these junctions remain to be understood.
Recent progress in scanning tunneling microscopy (STM) and spectroscopy (STS) technologies opens up a way to directly visualize the spatial variation of the electronic states in superconducting hetero-structures [3][4][5]. However, STM/STS studies on superconducting junctions made of unconventional superconductors are still demanding. There are two reasons which make it difficult to study unconventional junctions. First, it is often challenging to artificially fabricate well-defined junctions of unconventional superconductors. Second, in most of unconventional superconductors, surfaces are not electronically neutral; the resultant charge accumulation at the surfaces prevents STM/STS from accessing bulk superconducting properties. In this study, we solve these problems by inspecting the twin boundaries (TBs) in the nodal iron-based superconductor FeSe [6,7].
The TB is a crystallographic plane in a crystal shared by two neighboring domains with one being the mirror image of the other. The TBs are often formed by a tetragonal-to-orthorhombic structural phase transition, which reduces the four-fold (C 4 ) symmetry at high temperature to two-fold (C 2 ) symmetry at low temperature. In such a case, the orthorhombic crystal may contain the TBs parallel to the (110) plane, which act as an atomically well-defined junction. Some unconventional-superconductor-related materials, such as YBa 2 Cu 3 O 7−δ , AE(Fe 1−x Co x ) 2 As 2 (AE: alkaliearth element) and NaFeAs, do form TBs upon the tetragonal-to-orthorhombic transition which were identified by STM/STS measurements [8][9][10]. However, unavoidable surface state formation and/or insufficient amount of chemical doping prevent the STM/STS measurements to access superconductivity near TBs in these materials.
FeSe (superconducting transition temperature T c ≈ 9 K [11]) is a promising candidate for studying the ef-fects of TBs on unconventional superconductivity by STM/STS. Among various iron-based superconductors, FeSe has the simplest crystal structure [ Fig. 1(a)] in which electronically neutral two-dimensional FeSe layers are stacked along the c axis [11]. This guarantees the perfect cleaved surface which is electronically neutral. The tetragonal-to-orthorhombic structural phase transition, which is likely caused by the orbital ordering [12][13][14][15][16][17], occurs at T s ≈ 90 K and the TBs are spontaneously formed in the orthorhombic phase as illustrated in Fig. 1(b).
Band-structure calculations show that the Fermi surface of FeSe consists of hole cylinders around the zone center and compensating electron cylinders around the zone corner [18,19]. Several measurements, including penetration depth, quasiparticle interference, thermoelectric response [7], quantum oscillations [14,20,21], and angle-resolved photoemission spectroscopy (ARPES) [12][13][14][15]22] reveal that the Fermi surface in the orthorhombic phase consists of one hole and one (or two [14,21]) electron bands, both of which have very low carrier densities. The tunneling spectrum [6], temperature dependence of the penetration depth down to 80 mK and the residual thermal conductivity at T → 0 [7], all provide strong evidence that FeSe is an unconventional superconductor with line nodes in the superconducting gap.
The TBs in FeSe have been studied by low-temperature (4.2 K) STM/STS in the films grown by molecular beam epitaxy and the suppression of superconductivity by the TBs has been reported [23]. We performed STM/STS measurements at much lower temperature (≈ 0.4 K) in vapor-grown single crystals to examine the details of the superconducting gap and quasiparticle excitations near the TBs.

II. EXPERIMENTAL METHOD
STM/STS experiments were conducted in a constantcurrent mode with a commercial ultra-high vacuum verylow temperature STM (UNISOKU, USM-1300) modified by ourselves [24]. The samples used in this study were high-quality bulk single crystals grown using the vapor transport method [26]. Superconducting transition temperature defined at zero resistance is about 9 K. These crystals are undoped and stoichiometric, enabling us to investigate uniform and clean TBs. Samples were cleaved in-situ at liquid N 2 temperature to prepare clean and flat (001) surfaces. Immediately after cleaving, the samples were transferred to the STM unit kept below 10 K. We used electrochemically-etched polycrystalline tungsten wires for the scanning tips which were cleaned and sharpened in-situ by field evaporation using field-ion microscopy. The tunneling conductance g(r, E) ≡ dI t /dV s (r, E) reflecting the local density of states (LDOS) at a position r and energy E, was acquired by standard lock-in technique. Here, I t and V s denote the tunneling current and the sample-bias volt-age, respectively.

III. RESULTS AND DISCUSSION
A. Imaging the twin boundary in FeSe Figure 1(c) depicts an STM image of the cleaved surface of an FeSe single crystal at temperature T = 1.5 K. The image demonstrates the extremely small concentration of defects, i.e. about one defect per 5000 Fe atoms in the (001) plane. There is a shallow "groove" running along the [110] direction of the Fe lattice across which the unidirectional feature around the point defect is rotated by π/2, indicating that the "groove" represents the TB. We also observed that the elongated vortex cores [6], which were imaged by mapping g(r, E = 0) in a magnetic field, are rotated by π/2 across the TB [ Fig. 1(d)]. What is intriguing is that the vortices trapped at the TB are not elongated along the TB, demonstrating that the critical current density across the TB is comparable to that of the bulk. The STM image of the TB at a lower bias voltage is not a "groove" but a "ridge" [Figs. 1(e)]. This suggests that the apparent corrugations near the TB are primarily associated with the electronic-state variations; the actual surface topography near the TB may be essentially flat. A magnified STM image near the TB is shown in Fig. 1(f). A regular square lattice of the top-most Se atoms is well maintained even in the close vicinity of the TB. These observations indicate that the TB in FeSe is an atomically sharp superconducting junction with minimal strain to the lattice. (A detailed argument regarding the absence of strain is given in Appendix A.)

B. Local density of states across the twin boundaries
We examined the LDOS evolution across the TB by taking a series of g(r, E) along the line indicated in Fig. 2(a). Here and in the following, we are interested only in the evolution of g(r, E) along the x axis running perpendicular to the TB leaving the y coordinate constant, hence g(x, E). Figure 2(b) shows an intensity plot of g(r, E). Individual spectra taken at representative points are depicted in Fig. 2(c). At the position far away from the TB (I), g(r, E) exhibits a superconducting gap with clear quasiparticle peaks at ≈ ±2.5 meV. In addition to this main feature, there is a shoulder outside of the main peaks (≈ ±3.5 meV), which may represent multiple superconducting gaps [7]. In contrast to the case of fully-gapped superconductors in which g(r, E) = 0 in an extended E region near E = 0, g(r, E) in FeSe approaches zero only for E → 0 and apparently V-shaped, indicating the presence of line nodes [6]. Even right at the TB (III), the residual LDOS at E = 0 is negligibly small, indicating that the TB hardly gives rise to a pair breaking effect. In the vicinity of the TB, the quasi- Note that the pattern is rotated by π/2 between the two domains. The set-up conditions for imaging were Vs = +95 mV and It = 10 pA. (d) Zero-bias conductance image g(r, E = 0) at 1.5 K showing vortices. A magnetic field of 1 T was applied along the c axis. The tip was stabilized at Vs = +10 mV and It = 100 pA. A bias modulation amplitude V mod = 0.21 mVrms was used for spectroscopy. (e) A low-bias STM image at 1.5 K taken with Vs = +20 mV and It = 10 pA. The field of view for (c)-(e) is the same. (f) An atomic-resolution STM image near the TB which is running vertically in the center of the field of view. Vs = +95 mV and It = 100 pA.
particle peak and the shoulder associated with the superconducting gap diminish, and instead, sharp particlehole symmetric peaks appear at E ≈ ±1.5 meV. In the crossover region (II), the 1.5 meV peak coexists with the 2.5 meV peak, meaning that the former is not associated with the suppressed superconducting gap. The 1.5 meV peak diminishes within a distance of about 5 nm from the TB, which is close to the "averaged" in-plane coherence length ξ ab ≈ 5 nm obtained from the upper critical field H c2 ( c) ≈15 T [7,20]. These results suggest that the 1.5 meV peak represents the bound state induced by the TB.
Another interesting observation is that low-energy quasiparticle excitations are modified over a very long distance from the TB. High-resolution g(r, E) spectra at the positions of (I), (II) and (III) are plotted in Fig. 2

(d).
While the overall V-shaped behavior is maintained, the exact shape near the bottom of the gap depends on the position. In order to examine this behavior, we fit an empirical power-law g(r, E) ∝ |E| α to the low-energy (|E| < 0.5 meV) spectra and plot the exponent α as a function of the distance from the TB at x = 0 [ Fig. 2(e)]. Except close to the TB (|x| ξ ab ) where the 1.5 meV peaks dominate, α increases gradually with decreasing x by about ≈ 40%. This implies the suppression of the lowenergy quasiparticle excitations, most probably due to the opening of a small gap induced by the TB. The salient feature is that α continues to change even at |x| > 10ξ ab (≈ 50 nm), indicating an unexpectedly long-distance influence of the TB.
The long-distance TB effect on the LDOS can be seen in a more dramatic way in two junctions in series formed by two TBs. As shown in Fig. 3(a) we find an area where two TBs are running parallel to each other. The distance between the TBs is 34 nm, which is about 7 times larger than ξ ab . Figure 3(b) shows the spatial evolution of g(r, E) across the double TBs. Individual spectra at representative points are plotted in Fig. 3(c). The overall spectral features, the 2.5 meV peak, the 3.5 meV shoulder and the 1.5 meV peak observed near a single TB are all reproduced (positions I, II, and III). However, the lowenergy spectrum taken inside the central domain (position IV) shows a striking anomaly which is absent in the case of a single TB. Figure 3(d) depicts g(r, E) spectra at low energies. It is clear that, in between the double TBs, there is a finite energy range where g(r, E) is almost completely zero. The noticeable difference of the gap structure between inside and outside the central domain is clearly seen in Fig. 3(e), which shows the exponent α plotted as a function of the distance from one of the TBs; α is strongly enhanced in the central domain and peaks at the middle of the domain. The large power α ≈ 4, which is ≈ 3 times larger than the values at large x, essentially indistinguishable from an exponential energy dependence. This apparent large power again corroborates the finite gap opening in the excitation spectrum of quasiparticle.

C. Possible time-reversal-symmetry-broken state near the twin boundary
The above observations, the TB-induced bound states at finite energies and the suppression of the low-energy quasiparticle excitations over a length scale much longer than ξ ab , suggest a novel role of the TB in an unconventional superconductor. Before discussing the origin of these anomalies, we briefly review what can be expected at a TB of FeSe. Recent high-resolution laser-ARPES measurements of FeSe indicate that the hole cylinder is fully gapped [27], implying that the line nodes are present on the electron cylinder, as has been also inferred from vortex imaging [6]. Given this information, we consider two possible phase structures for symmetry of the superconducting gap across a TB as illustrated in Fig. 4, where either the global phase of the superconducting gap is fixed to the crystallographic axis [ Fig. 4(a)] or is flipped across the TB [ Fig. 4(b)]. It should be noted that the sign of either the nodal gap or the nodeless gap is reversed between the two domains in Fig. 4(a) or Fig. 4(b), respectively. This means that the amplitude of at least one of the gaps vanishes at the TB, giving rise to the zero-energy quasiparticle state that should appear as a zero-energy peak in g(r, E). This argument applies not only for the particular phase structure shown in Fig. 4 but also for a general case in which nodal and nodeless gaps reside on multiple Fermi surfaces.
The observed bound-state peak at 1.5 meV apparently contradicts this conjecture and suggests instead that the TB induces an additional gap component which shifts the position of a zero-energy peak to a finite energy. We point out that, as long as the induced gap is real, a sum of the bulk gap and the TB-induced gap reverses its sign at a finite distance from the TB and still gives rise to a zero-energy peak. However, as shown in Fig. 2(b), we did not observe a zero-energy peak in g(r, E) over more than 100 nm from the TB. Thus, we speculate that the induced gap has an imaginary component, which means that time reversal symmetry is locally broken near the TB. In such a case, bound state peaks are located at finite energies E = ±∆ cos(δϕ/2) because the phase shift δϕ on the TB is reduced from π [28,29]. Here, ∆ is the amplitude of the superconducting gap. The possibility of the TB-induced time-reversal-symmetry-broken state has been argued theoretically in d-wave YBa 2 Cu 3 O 7−δ with a small s-wave component [30], but the experimental observation is still lacking.
In order to substantiate the relevance of this scenario, we have calculated the spatial evolution of the LDOS for a model order parameter with broken time-reversal symmetry near TBs. The C 2 -symmetric order parameter is represented by a sum of the isotropic component ∆ iso and the four-fold nodal component ∆ 4φ sin(2φ), where φ is the azimuthal angle in the momentum space; see Fig. 4. We assume that the global phase of the order parameter is fixed to the crystallographic axis, that is, the nodal component changes its sign across a TB as shown in Fig. 4(a). The spatial variation of ∆ 4φ around a TB at x = x 0 is modeled by the form where the x axis is taken to be perpendicular to the TB, and ∆ bulk 4φ is the amplitude of ∆ 4φ in the bulk. The phase The phase variable θ(x) is assumed to take a nonvanishing value near the TB and exponentially decay with another length scaleξ. It is important to note that the characteristic lengthξ for the local timereversal symmetry breaking can be much longer than the coherence length ξ [30]. (The derivation of the length ξ is given in Appendix B.) To account for low-energy excitations near the nodes observed in the LDOS, we focus on the electron cylinder with nodal gaps by setting the parameters ∆ iso = 0.2∆ 0 and ∆ bulk 4φ = 0.8∆ 0 . As a model order parameter with a TB at x = x 0 = 0, we take θ(x) = (π/6)sech(x/ξ) withξ = 5ξ, which gives ∆ 4φ (x = 0) = (i/2)∆ bulk 4φ . The order parameter ∆ 4φ (x) is shown in Fig. 5(a), where |∆ 4φ | changes with the length scale ξ while Im(∆ 4φ ) decays with the longer length scale FIG. 4: Schematic illustration of the phases of the superconducting gaps across the TB shown by the red line. Top panel represents the iron lattice near the TB together with the momentum-space phase structure of the superconducting gaps opening at multiple Fermi cylinders, a hole cylinder at the center and electron cylinders at the corner of the Brillouin zone (black broken square). Different colors (red and blue) denote different signs of the phase. We assume that the gap node exists on the electron cylinder and the sign reversal is between the main lobe of the gap on the electron cylinder and the gap on the hole cylinders but the argument given in the text applies not only for this particular case but also for other cases. There are two possibilities: either the phase structure is fixed to the lattice (a) or is flipped across the TB (b). In the former case, the nodal component ∆ 4φ should change its sign across the TB, whereas the sign of the isotropic component ∆iso (either due to the fully gapped Fermi cylinder or associated with the C2 symmetry of the nodal gap) should be reversed in the latter case.
ξ. The phase ϕ abruptly changes near the TB and gradually approaches 0 or π. Using this order parameter, we calculate the spatial dependence of the LDOS within the quasi-classical approximation [31]. Figure 5(b) shows the global peak structure of the LDOS at representative points calculated with energy smearing of η = 0.03∆ 0 . Far from the TB, namely in the bulk (I), the LDOS has peaks at |E| = ∆ bulk 4φ ± ∆ iso . On the TB (III), the peaks observed in the bulk are suppressed, and alternative peaks appear at E ≈ ±0.4∆ 0 , which correspond to the bound states whose energies are shifted from E = 0 due to the local time-reversal symmetry breaking in ∆ 4φ . The bound-state peaks disappear at x = 3ξ (II) since their wave functions decay into the bulk with the length scale ξ. These features of the calculated LDOS arising from the bound states are consistent with the LDOS peaks observed at E ≈ ±1.5 meV by STM. We show in Fig. 5(c) the LDOS at lower energy scale which has been calculated with a much smaller smearing factor η = 0.001∆ 0 . The clear V-shaped LDOS in the bulk (I) changes to the U-shaped LDOS upon approaching the TB, in agreement with the increase of the exponent α evaluated from the experimental data [ Fig. 2(e)]. The low-energy LDOS is finite at x = 0 (III) and x = 3ξ (II) because low-energy quasiparticles with momenta along the nodal directions of the bulk gap can linger over long distance and reach the TB, even though the local gap, does not vanish near the TB where Im(∆ 4φ ) = 0. We also calculate the LDOS for double TBs located at x = ±x 0 = ±3.5ξ, taking the model order parameter of the form shown in Fig. 5(d). We assume that the distance 2x 0 between the TBs is in the range ξ ≪ x 0 ξ . The phase θ(x) is an even function of x and takes a maximum value at x = 0. The global peak structure of the LDOS and its low-energy blowup are shown for representative points along the x direction in Fig. 5(e) and 5(f), respectively. The large peaks at |E| ≈ 0.4∆ 0 on a TB (III) and the small peaks at |E| ≈ 0.2∆ 0 at the middle point x = 0 between the TBs (IV) in Fig. 5(e) originate from the same dispersive mode of bound states at a TB. The calculated LDOS spectrum at x = 0 between the double TBs (IV) in Fig. 5(f) exhibits a clear energy gap extending over the region |E| 0.1∆ 0 , reflecting the existence of a larger local gap at x = 0, where the bulk low-energy quasiparticles cannot reach. We conclude that the local gap enhanced by the local time-reversal symmetry breaking near TBs over the length scaleξ can explain the strong suppression of the LDOS between the two TBs observed in our STM/STS experiments.

IV. SUMMARY
We have reported on the visualization of the atomic scale variation of the quasiparticle states of the nodal superconductor FeSe near TBs that enforce a sign inversion at least one of the superconducting gaps opening on multiple Fermi cylinders. In contrast to the expectation that the sign inversion generates a zero-energy quasiparticle bound state near the TB, the TB-induced quasiparticle states are not at zero but at finite energies E ≈ ±1.5 meV. Moreover, the low-energy excitation spectrum is affected by the TB over an extremely long distance, which is a few tens of times larger than ξ ab . An even more dramatic change in the low-energy spectrum has been detected in the region between double TBs separated by a distance ≈ 7ξ ab , where the quasiparticle weight near the Fermi energy is almost completely removed in the energy range |E| 0.2 meV. These observations are qualitatively reproduced by a phenomenological model which assumes that the TB induces locally a superconducting state that breaks time-reversal symmetry.
Our results suggest several important directions for future studies. A microscopic mechanism that induces the time-reversal symmetry broken state is elusive and should be investigated theoretically. It is also interesting to go beyond the quasi-classical approximation because FeSe is a unique superconductor whose Fermi energy is of the same order as the superconducting gap, placing this system in the BCS-BEC crossover regime [7]. Experiments that directly probe the time-reversal symmetry breaking, such as muon-spin rotation and local magnetometry, are highly desirable and would give us further insights into the unconventional superconducting junctions. We anticipate that the TBs in FeSe will stimulate further research on the role of the phase of the superconducting order parameter near the interface, which has been difficult to access experimentally.

Acknowledgments
This work has been supported by Japan -Germany Research Cooperative Program, KAKENHI from JSPS and Project No. 56393598 from DAAD, and the "Topological Quantum Phenomena" (No. 25103713) KAK-ENHI on Innovative Areas from MEXT of Japan.

Appendix A: Absence of lattice distortion induced by the twin boundary
Although STM has a high spatial resolution, possible creep in the piezoelectric scanner and/or the thermal drift make it difficult to estimate the small distortions in the topographic image. Here we utilize the so-called Lawler-Fujita algorithm [32] to deduce the lattice distortion and show that the TB-induced strain is negligibly small.
First we briefly explain the principle of the methodology. The observed STM topography T (r), which mainly represents the top-most Se lattice, can be expressed as Here, T 0 is the amplitude of the Se-lattice modulation, q x and q y are wave vectors for the Se lattice, and · · · represent all other modulations. The distortions from the perfect lattice is described by the displacement field u(r) that can be regarded as a spatially varying phase of the q x and q y modulations. This approximation is justified as long as the length scale of distortions is much longer than the Se-Se distance a Se . Standard phase-sensitive detection scheme can be used to evaluate u(r). By multiplying T (r) and the reference signal cos (q x · r), we get T (r) cos (q x · r) = T 0 2 cos (q x · u(r)) + cos (2q x · r − q x · u(r)) + cos ((q x + q y ) · r − q y · u(r)) + cos ((−q x + q y ) · r − q y · u(r)) All terms except the first exhibit periodic spatial modulations, which can be removed by low-pass Fourier filtering LPF {· · · }.
By using the quadrature reference sin (q x · r), we get Therefore, we obtain u x (r), the x component of u(r) as The y component u y (r) can also be deduced as u y (r) = a Se 2π tan −1 LPF {T (r) sin (q y · r)} LPF {T (r) cos (q y · r)} . (A6) A schematic model of atomic arrangement near the TB is shown in Fig. 6(a). We expect that the orthorhombic distortion affects the atomic arrangement along the y direction across the TB, while the periodicity along the x direction remains intact. In order to verify this model and to check if there is an additional lattice distortion, we calculated u x (r) and u y (r) of the high-resolution STM image containing a TB running along the y direction [Fig 6(b)]. Reference wave vectors q x and q y were obtained by the Fourier analysis in the left domain. For lowpass Fourier filtering, we picked up only long-wavelength components by using a Gaussian mask with half width at the half maximum of 0.21(2π/a Se ). Since there is a translational symmetry along the TB, we average u x (r) and u y (r) along the y direction, yielding u avg x (x) and u avg y (x), respectively. This significantly enhances the signal-tonoise ratio.
Figures 6(c) and (d) show u avg x (x) and its x derivative. There is no noticeable anomaly in both u avg x (x) and du avg x (x)/dx, except for the smooth background associated with the creep of the piezoelectric scanner. By contrast, u avg y (x) exhibits a sharp kink at the TB [ Fig. 6(e)]. These features are consistent with the model shown in Fig. 6(a). It should be noted that du avg y (x)/dx shown in Fig. 6(f) is almost completely constant in both domains, indicating that the TB-induced strain to the lattice is negligibly small. The observed value of du avg y (x)/dx ≈ −1.1 × 10 −2 in the right domain means that the angle β defined in Fig. 6(a) is +0.63 • . This means that orthorhombic distortion (b − a)/(b + a) ≈ 2.8 × 10 −3 , being consistent with the X-ray result [33]. Even if there were an additional lattice distortion associated with the TB, it should be much smaller than this tiny orthorhombic distortion which we have clearly detected.
Appendix B: Asymptotic forms of the order parameter derived by the Ginzburg-Landau theory We derive asymptotic forms of the order parameter far from TBs using the Ginzburg-Landau (GL) theory. We consider the GL free-energy functional for tetragonal symmetric systems [30] as an expansion in the isotropic s-wave component ∆ iso and the four-fold d-wave component ∆ 4φ of the order parameter: where we have neglected the vector potential as it does not play an important role in our discussion. The coefficients b µ , K µ , and K are positive andã µ (T ) = a µ (T /T cµ − 1) with positive a µ . The differential operator ∇ = (∂ a , ∂ b ) is defined according to the crystal axes a and b. As in Ref. 30, we assume γ 2 > 0, so that the free energy is minimized at ϕ = ±π/2, and the timereversal-symmetry-broken s ± id state is stabilized when both ∆ iso and ∆ 4φ are finite. The effect of orthorhombic distortion is taken into account by adding the following term to the free-energy functional [30]: where c is a positive parameter and ǫ = ǫ aa − ǫ bb is the parameter of the orthorhombic lattice distortion. The total free energy for a uniform state in the bulk is then given by where ∆ µ = |∆ µ |e iϕµ and the relative phase ϕ = ϕ 4φ − ϕ iso . If c|ǫ| ≥ 2γ 2 |∆ iso ||∆ 4φ |, then the free energy is minimized at ϕ = 0 for ǫ < 0 and at ϕ = π for ǫ > 0. In the following discussion we assume that this inequality is satisfied and the time reversal symmetric s ± d state is realized in the bulk. Next, we consider a TB located at x = x 0 along the y axis, where the x and y axes are rotated by 45 • from the crystalline axes, x = (a − b)/ √ 2 and y = (a + b)/ √ 2. The orthorhombic lattice distortion parameter ǫ changes its sign across the TB. We assume ǫ(x) → ∓|ǫ| for x → ±∞, so that the s ± d state is realized in x → ±∞. Near the TB where ǫ(x) is small, the s ± id state is favored. Then, the area density of the total free energy is given by Let us first assume that ∆ iso is a real and uniform order parameter while ∆ 4φ changes its sign across the TB, as shown in Fig. 4(a). If we restrict ∆ 4φ to be real, then ∆ 4φ varies over the coherence length [30] ξ = K 4φ a 4φ + 6b 4φ |∆ bulk 4φ | 2 + (γ 1 + γ 2 )|∆ iso | 2 , where |∆ bulk 4φ | is the amplitude of ∆ 4φ in the bulk. However, we expect that time-reversal symmetry should be locally broken at the TB. Thus, we allow ∆ 4φ to be complex, ∆ 4φ (x) = |∆ 4φ (x)|e iϕ(x) . With this order parameter, the total free energy is given by In the bulk region (x − x 0 ≫ ξ) where ∂ x |∆ 4φ | = 0 and ǫ(x) = −|ǫ|, the GL differential equation to minimize Since ϕ ≪ 1 far away from the TB, we can linearize the differential equation and find the relative phase to decay as ϕ ∝ exp −x/ξ with the characteristic length ξ = K 4φ |∆ bulk 4φ | |∆ iso |(c|ǫ| − 2γ 2 |∆ iso ||∆ bulk 4φ |) .