**1. Introduction**

The discovery of iron-based materials with a high superconducting transition temperature has triggered great interest for both fundamental studies and practical applications in this field [1]. Among these superconductors, the FeSe system has a simple crystal structure, clean superconducting phase and low toxicity, making it an appealing candidate for studying the superconducting properties of Fe-based compounds [2]. The FeSe layers consist of square lattices of Fe atoms with tetrahedrally coordinated covalent bonds to the Se anions, and the lattice constant perpendicular to the layered plane is about 0.55 nm [3]. The Fermi surface of this compound consists of one electron and one-hole thin cylinders through Shubnikov–de Haas oscillations [4]. At around 100 K, FeSe shows a structural phase transition from tetragonal to orthorhombic without an accompanying magnetic phase transition and becomes superconducting below 8 K [2,5–7]. The picture of two-gap superconductivity has been clearly confirmed by the scanning tunneling microscopy measurements, multiple Andreev reflection spectroscopy, specific heat measurements and other experiments [8–12].

High-quality, superconducting thin films have an important role in applications and basic research of superconductivity. In this respect, preparing high-quality thin-film samples not only satisfies the demands for some measurements of basic physical properties but also provides suitable bases for making tunneling junctions, which determines several important superconducting parameters, such as gap value and paring symmetry [13]. Previous studies examining the thickness dependence of FeSe have been limited to measurements on thin films grown using techniques such as molecular beam expitaxy [14,15], pulsed laser deposition [16] and radio-frequency sputtering [17], all of which require welloptimised growth protocols. An alternative to the growth of thin films is to create devices by mechanical exfoliation of high-quality single crystals. With this method, a series of FeSe superconducting films with thicknesses ranging from 470 to 2.2 nm was obtained on the Si/SiO2 substrate [18]. A dramatic depression of *Tc* has been observed when the thickness is smaller than 27 nm. Farrar et al. have also observed similar behavior: a sharp decrease in superconductivity occurs at *d* < 25 nm [19]. One possible origin for this *Tc* suppression could be the increase in disorder scattering with reducing thickness [19,20]. Another alternative possibility is due to the interaction between FeSe thin film and the

**Citation:** Ye, C.; Che, J.; Huang, H. Boundary Effect and Critical Temperature of Two-Band Superconducting FeSe Films. *Crystals* **2023**, *13*, 18. https://doi.org/ 10.3390/cryst13010018

Academic Editor: Jacek Cwik ´

Received: 2 December 2022 Revised: 18 December 2022 Accepted: 18 December 2022 Published: 22 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

substrate [21]. However, up to now, there is still no consensus on the explanation of the experimental data mentioned above.

In this paper, we propose that the *Tc* dependence on the film thickness is due to the influence of the boundary effect between the two-band superconductor and the insulator (or vacuum). We first introduce the appropriate boundary conditions for the Ginzburg– Landau (GL) order parameters at the superconductor–insulator interface. We then give a microscopic analysis of these new boundary terms based on two-band Bogoliubov–de Gennes theory. Our theoretical result is consistent with the experimental data of FeSe films, which suggests the boundary effect is an important factor for the understanding of superconducting properties in this iron-based compound.

The rest of this article is structured as follows. In Section 2, we first review two-band Bogoliubov–de Gennes theory and GL equations, and then give a microscopic derivation of proper boundary conditions for the superconductor-insulator interface. In Section 3, we perform the calculations on the film-thickness-dependence of critical temperature for the FeSe compound in the context of GL theory. Finally, we conclude this article in Section 4.

#### **2. Theoretical Scheme**

Based on the previous literature [22–31], we can write the Hamiltonian of a two-band superconductor as

$$H = \sum\_{i\sigma} c\_{i\sigma}^{\dagger}(r)\hat{h}(r)c\_{i\sigma}(r) - \sum\_{i i'} g\_{i i'} c\_{i\uparrow}^{\dagger}(r)c\_{i\downarrow}^{\dagger}(r)c\_{i'\downarrow}(r)c\_{i'\uparrow}(r),\tag{1}$$

where *i*, *i* <sup>=</sup> 1, 2 are the band indices and *<sup>σ</sup>* <sup>=</sup>↑, <sup>↓</sup> is the spin index. <sup>ˆ</sup> *<sup>h</sup>*(*r*) is the single-particle Hamiltonian of the normal metal, and *gii* are the electron–phonon interaction constants with *g*<sup>12</sup> = *g*21.

We can introduce the gap functions as

$$\Psi\_i(r) = -\sum\_{i'} g\_{ii'} \langle c\_{i'\downarrow}(r) c\_{i'\uparrow}(r) \rangle \tag{2}$$

and transform the Hamiltonian into the mean-field form

$$H\_{eff} = \sum\_{i\sigma} c\_{i\sigma}^{\dagger}(r)\hat{h}(r)c\_{i\sigma}(r) + \sum\_{i} [\Psi\_i(r)c\_{i\uparrow}^{\dagger}(r)c\_{i\downarrow}^{\dagger}(r) + \text{H.c.}].\tag{3}$$

This effective Hamiltonian can be diagonalized by means of the Bogoliubov transformation with *b* and *b*† as the annihilation and creation operators of quasi-particle excitations:

$$c\_{i\uparrow}(r) = \sum\_{k} [\mu\_{ik}(r)b\_{ik\uparrow} - v\_{ik}^\*(r)b\_{ik\downarrow}^\dagger] \tag{4}$$

and

$$c\_{i\downarrow}(r) = \sum\_{\mathbf{k}} [u\_{i\mathbf{k}}(r)b\_{i\mathbf{k}\downarrow} + v\_{i\mathbf{k}}^{\*}(r)b\_{i\mathbf{k}\uparrow}^{\dagger}] \tag{5}$$

where *k* is the wave vector. As a result, the effective Hamiltonian can be written as

$$H\_{eff} = E\_{\mathfrak{F}} + \sum\_{ikr} E\_{ik} b\_{ikr}^{\dagger} b\_{ikr} \,\tag{6}$$

*Eg* is the ground state energy, and *Eik* is the energy of the excitation.

Using the commutator [*ciσ*(*r*), *Heff*], together with Equations (4)–(6), we can obtain the Bogoliubov–de Gennes equations for a two-band superconductor [32–35]

$$
\begin{pmatrix}
\hat{h} & \Psi\_i(\mathbf{r}) \\
\Psi\_i^\*(\mathbf{r}) & -\hat{h}^\*
\end{pmatrix}
\begin{pmatrix}
u\_{ik}(\mathbf{r}) \\
v\_{ik}(\mathbf{r})
\end{pmatrix} = E\_{i\mathbf{k}} \begin{pmatrix}
u\_{ik}(\mathbf{r}) \\
v\_{ik}(\mathbf{r})
\end{pmatrix}.
\tag{7}
$$

From Equation (6), we can also obtain *b*† *<sup>i</sup>k*<sup>↑</sup>*bik*<sup>↑</sup> <sup>=</sup> *<sup>f</sup>*(*Eik*) with *<sup>f</sup>*(*Eik*)=[<sup>1</sup> <sup>+</sup> exp(*Eik*/*kBT*)]−1. Then, with Equation (2), we can transform the self-consistent gap equations into

$$\Psi\_i(\mathbf{r}) = \sum\_{i'k} g\_{ii'} v\_{i'k}^\*(\mathbf{r}) u\_{i'k}(\mathbf{r}) \times [1 - 2f(E\_{i'k})].\tag{8}$$

In analogy with the single-band case [36], for small gap functions Ψ*i*, we can obtain the linearized form of self-consistency conditions from Equations (7) and (8) as

$$\Psi\_i(\mathbf{r}) = \sum\_{i'} \int \mathcal{K}\_{ii'}(\mathbf{r}, \mathbf{r'}) \Psi\_{i'}(\mathbf{r'}) d\mathbf{r'} \tag{9}$$

with the kernel

$$K\_{ii'}(r,r') = \mathcal{g}\_{ii'}k\_BT \sum\_{\mathbf{k}\mathbf{k'}} \sum\_{\omega} \frac{\Phi\_{i'\mathbf{k}}^{\*}(r')\Phi\_{i'\mathbf{k'}}^{\*}(r')\Phi\_{i'\mathbf{k}}(r)\Phi\_{i'\mathbf{k'}}(r)}{(\varepsilon\_{i'\mathbf{k}} - i\hbar\omega)(\varepsilon\_{i'\mathbf{k'}} + i\hbar\omega)}.\tag{10}$$

<sup>Φ</sup>*ik*(*r*) is defined as the normal-state eigenfunction of the electron; <sup>ˆ</sup> *<sup>h</sup>*Φ*ik* <sup>=</sup> *<sup>ε</sup>ik*Φ*ik*. The frequency *ω* = (2*ν* + 1)*πkBT*/¯*h*, where *ν* is an integer.

With the explicit expressions of the kernels in the bulk system and the addition of nonlinear terms to the gap equations, we can obtain the two-band GL equations from Equation (9) as [27]

$$
\alpha\_1(T)\Psi\_1 + \beta\_1|\Psi\_1|^2\Psi\_1 - \gamma\_1\nabla^2\Psi\_1 - \epsilon\_{12}\Psi\_2 = 0\tag{11}
$$

and

$$
\alpha\_2(T)\Psi\_2 + \beta\_2|\Psi\_2|^2\Psi\_2 - \gamma\_2\nabla^2\Psi\_2 - \epsilon\_{12}\Psi\_1 = 0\tag{12}
$$

with the GL parameters

$$\alpha\_{1,2} = N\_{1,2} \left[ \frac{\lambda\_{22,11}}{\lambda} - \frac{1}{\lambda\_{\text{max}}} - \ln \left( \frac{T\_{c0}}{T} \right) \right], \quad \beta\_i = \frac{7\zeta(3)N\_i}{16\pi^2 (k\_B T\_{c0})^2} \tag{13}$$

$$\gamma\_i = \frac{7\zeta(3)\hbar^2 N\_i v\_{Fi}^2}{16\pi^2 (k\_B T\_{c0})^2} \text{ and } \epsilon\_{12} = \frac{N\_1 \lambda\_{12}}{\lambda} = \frac{N\_2 \lambda\_{21}}{\lambda}. \tag{14}$$

*λii* = *giiNi* with *Ni* being the density of states at the Fermi level for each band; *λ* = *<sup>λ</sup>*11*λ*<sup>22</sup> <sup>−</sup> *<sup>λ</sup>*12*λ*<sup>21</sup> and *<sup>λ</sup>*max <sup>=</sup> <sup>1</sup> 2 (*λ*<sup>11</sup> <sup>+</sup> *<sup>λ</sup>*22) + (*λ*<sup>11</sup> <sup>−</sup> *<sup>λ</sup>*22)<sup>2</sup> <sup>+</sup> <sup>4</sup>*λ*12*λ*21 are the determinant and the largest eigenvalue of *λ*-matrix, respectively. *Tc*<sup>0</sup> is the bulk critical temperature, and *vFi* is the average Fermi velocity for each band.

In the spatially homogeneous case, we can neglect the gradient *γ*-terms. Equations (11) and (12) will yield the gap equations at *T* = *Tc*0:

$$
\lambda \begin{pmatrix} \lambda\_{11} & \lambda\_{12} \\ \lambda\_{21} & \lambda\_{22} \end{pmatrix} \begin{pmatrix} \Psi\_1 \\ \Psi\_2 \end{pmatrix} = \lambda\_{\text{max}} \begin{pmatrix} \Psi\_1 \\ \Psi\_2 \end{pmatrix} . \tag{15}
$$

which obviously give a consistent result.

Meanwhile, we can write down the boundary conditions for the two-band GL theory at the interface between the two-band superconductor and insulator (or vacuum) as

$$\nabla \Psi\_i \cdot \mathbf{s}|\_S = -\sum\_{i'} A\_{ii'} \Psi\_{i'} \tag{16}$$

with *S*, the boundary coordinate, and *Aii* , as some constants. From Equation (16), we can see that the ordinary Neumann boundary condition corresponds to *Aii* = 0. However, we will show that the boundary effect induced by these *A*-terms is important for the understanding of the *Tc* suppression in FeSe thin films.

At this stage, we would also like to point out the boundary effect and different interband interactions in multi-band superconductors have already been extensively studied in the literature. Babaev et al. presented a microscopic study on the behavior of the order parameters for two-band superconductors based on the free boundary condition [37]. With the Neumann boundary condition, the stable edge states and the dynamic response of such states to an external applied current have been investigated in the timedependent Ginzburg–Landau formalism for the two-band mesoscopic superconductors [38]. Aguirre et al. also discussed the effects of different interband interactions on the vortex states by solving the two-band Ginzburg–Landau equations with the Neumann boundary condition [39,40]. However, to explain the suppression of critical temperature with the decrease in film thickness, we need to conduct detailed microscopic analysis and derive the correct boundary terms based on the two-band Bogoliubov–de Gennes theory for the superconductor FeSe.

Now, we try to derive these boundary *A*-terms in Equation (16) from the two-band Bogoliubov–de Gennes theory. In all cases, we assume that there is no current flowing through the boundary. The equation to be solved reads

$$\Psi\_i(\mathbf{s}) = \sum\_{\vec{i}'} \int \mathcal{K}\_{\vec{i}\vec{\ell}'}(\mathbf{s}, \mathbf{s}') \Psi\_{\vec{i}'}(\mathbf{s}') d\mathbf{s}' \tag{17}$$

where *s* measures the normal distance from the boundary. For simplicity, we set the crosssection of the boundary as 1. *Kii*(*s*,*s* ) is defined by Equation (10), and due to the existence of interface Ψ*i*(*s*) will decrease exponentially in the insulating regime.

Following the procedure suggested by de Gennes [36], we suppose that the form of gap functions close to the surface behaves as

$$\Psi\_i(s) = \Psi\_{i0} + \left(\sum\_{i'} A\_{ii'} \Psi\_{i'0}\right)s \tag{18}$$

with Ψ*i*<sup>0</sup> being the gap function at the boundary and *s* > 0 inside the superconductor. It is easy to see that the boundary condition in Equation (16) follows naturally from Equation (18). However, beyond the scale of the coherence length from the boundary, the linear dependence definitely becomes invalid. Ψ*<sup>i</sup>* will then have a negative curvature and reach the BCS value deep in the superconductor.

If we introduce *K*<sup>0</sup> *ii*(*s*,*s* ) as the kernel of gap functions in the bulk metal, we can then transform Equation (17) as

$$\begin{split} \Psi\_{i}(\mathbf{s}) - \sum\_{i'} \int K\_{ii'}^{0}(s, s') \Psi\_{i'}(s') ds' \\ = - \sum\_{i'} \int [K\_{ii'}^{0}(s, s') - K\_{ii'}(s, s')] \Psi\_{i'}(s') ds' \equiv - \sum\_{i'} H\_{ii'}(\mathbf{s}). \end{split} \tag{19}$$

From Equations (11) and (12) with the higher order *β*-terms omitted, while also noting that *K*<sup>0</sup> *ii*(*s*,*s* ) = *K*<sup>0</sup> *ii*(*s* − *s* ) due to the translational symmetry, we can read out the Laplace transformation of *K*<sup>0</sup> *ii* close to the critical temperature as

$$K\_{ii'}^{0}(p) = \frac{\lambda\_{ii'}}{\lambda\_{\text{max}}} + \frac{\lambda\_{ii'}\gamma\_{i'}}{N\_{i'}}p^2. \tag{20}$$

By plugging Equation (20) into Equation (19), we can get

$$\Psi\_i(p) - \sum\_{\vec{i}'} (\lambda\_{i\vec{i}'} / \lambda\_{\text{max}}) \Psi\_{\vec{i}'}(p) - \sum\_{\vec{i}'} (\lambda\_{i\vec{i}'} \gamma\_{\vec{i}'} / N\_{\vec{i}'}) p^2 \Psi\_{\vec{i}'}(p) = -\sum\_{\vec{i}'} H\_{\vec{i}\vec{i}'}(p). \tag{21}$$

Ψ*i*(*p*) and *Hii*(*p*) are the Laplace transformations of Ψ*i*(*s*) and *Hii*(*s*), respectively. Since the first two terms of the left-handed side in Equation (21) can be approximately canceled out according to Equation (15), we have

$$\sum\_{i'} (\lambda\_{ii'} \gamma\_{i'} / N\_{i'}) p^2 \Psi\_{i'}(p) = \sum\_{i'} H\_{ii'}(p). \tag{22}$$

We can see that both sides in Equation (22) take the main contribution from the boundary region. Notice that the Laplace transformation of the gap functions in Equation (18) takes the form

$$
\Psi\_i(p) = \frac{\Psi\_{i0}}{p} + \sum\_{i'} \frac{A\_{ii'} \Psi\_{i'0}}{p^2}.\tag{23}
$$

Then, at *p* → 0, we will obtain from Equation (22)

$$\sum\_{i'i''} (\lambda\_{ii'} \gamma\_{i'} / N\_{i'}) A\_{i'i''} \Psi\_{i''0} = \sum\_{i'} H\_{ii'} (p=0). \tag{24}$$

According to de Gennes' analysis [36,41], from the sum rules

$$\int K\_{\rm ii'}^{0}(s,s')ds' = \frac{\lambda\_{\rm ii'}}{\lambda\_{\rm max}} \text{ and } \int K\_{\rm ii'}(s,s')ds' = \frac{\lambda\_{\rm ii'}N\_{\rm i'}(s)}{\lambda\_{\rm max}N\_{\rm i'}}\tag{25}$$

with *Ni*(*s*) as the local density of states at the Fermi surface, we can write the Laplace transformation of the kernel difference at *p* → 0 as

$$H\_{\rm i\'}(p=0) = \int H\_{\rm i\'}(s)ds = \frac{\lambda\_{\rm i\'}\Psi\_{\rm i'0}}{\lambda\_{\rm max}} \int \frac{\Psi\_{\rm i'}(s)}{\Psi\_{\rm i'0}} \left[1 - \frac{N\_{\rm i'}(s)}{N\_{\rm i'}}\right] ds. \tag{26}$$

Ψ*i*(*s*)/Ψ*i*<sup>0</sup> approaches zero in the insulating region, and is of the order of one in the metallic region. *Ni*(*s*)/*Ni* also passes from 0 → 1 a few atoms from the boundary. Therefore, the integrand in Equation (26) is nonvanishing only in a width of the order of the lattice constant *a*. We can then estimate *Hii*(*p* = 0) as

$$H\_{\text{ii'}}(p=0) = \frac{\lambda\_{\text{ii'}}a}{\lambda\_{\text{max}}} \Psi\_{i'0}.\tag{27}$$

By combining Equation (24) with Equation (27), we can finally obtain

$$A\_{ii} = \frac{N\_i a}{\gamma\_i \lambda\_{\text{max}}} \text{ and } A\_{12} = A\_{21} = 0. \tag{28}$$

With these formulae, we successfully demonstrate the microscopic origin of boundary conditions in Equation (16).

Based on two-band GL theory, we can write the supercurrent at the boundary *S* as

$$J\_S = -ie\hbar \sum\_{i=1}^{2} \frac{1}{m\_i} [\Psi\_i^\* (\nabla \Psi\_i \cdot \mathbf{s})|\_S - \Psi\_i (\nabla \Psi\_i^\* \cdot \mathbf{s})|\_S]. \tag{29}$$

According to our boundary conditions, we have the supercurrent

$$J\_S = -\frac{\mathrm{i}e\hbar}{m\_1}[\Psi\_1^\*(-A\_{11}\Psi\_1) - \Psi\_1(-A\_{11}\Psi\_1^\*)] - \frac{\mathrm{i}e\hbar}{m\_2}[\Psi\_2^\*(-A\_{22}\Psi\_2) - \Psi\_2(-A\_{22}\Psi\_2^\*)] = 0,\tag{30}$$

which definitely gives a consistent result.

#### **3. Critical Temperature of FeSe Films in Ginzburg–Landau Theory**

In this section, we try to understand the film-thickness-dependence of *Tc* for FeSe based on the boundary effect mentioned above. We suppose that the film extends from *z* = −*d*/2 to *z* = *d*/2 and the film thickness is *d*.

From Equations (11) and (12), two-band GL equations can be written as

$$
\begin{pmatrix}
\hat{H}\_{11} & \hat{H}\_{12} \\
\hat{H}\_{21} & \hat{H}\_{22}
\end{pmatrix}
\begin{pmatrix}
\Psi\_1(\mathbf{r}) \\
\Psi\_2(\mathbf{r})
\end{pmatrix} = \mathbf{0} \tag{31}
$$

with

$$
\hat{H}\_{\text{ii}} = -\gamma\_i \nabla^2 + a\_i(T) \tag{32}
$$

and

$$
\hat{H}\_{12} = \hat{H}\_{21} = -\epsilon\_{12}.\tag{33}
$$

Note that close to the critical temperature, the magnitudes of order parameters are small, and the higher order *β*-terms can be neglected.

Similarly to the single-band case [41], we set the form of gap functions for the superconducting film as

$$
\begin{pmatrix}
\Psi\_1(z) \\
\Psi\_2(z)
\end{pmatrix} = \begin{pmatrix}
\eta\_1 \cos(k\_1 z) \\
\eta\_2 \cos(k\_2 z)
\end{pmatrix} \tag{34}
$$

with *η<sup>i</sup>* as the constant. Then, from the boundary conditions

$$\left. \frac{d\Psi\_i}{dz} \right|\_{z=\pm d/2} = \mp A\_{ii} \Psi\_{i\nu} \tag{35}$$

we have *ki*, satisfying

$$k\_i \tan\left(\frac{k\_i d}{2}\right) = A\_{ii}.\tag{36}$$

Let us introduce

$$H\_{\vec{i}\vec{i}'} = \langle \Psi\_i | \hat{H}\_{\vec{i}\vec{i}'} | \Psi\_{\vec{i}'} \rangle = \int\_{-d/2}^{d/2} \Psi\_i(z) \hat{H}\_{\vec{i}\vec{i}'} \Psi\_{\vec{i}'}(z) dz. \tag{37}$$

We can transform Equation (31) into

$$
\left( \begin{array}{cc} H\_{11} & H\_{12} \\ H\_{21} & H\_{22} \end{array} \right) \left( \begin{array}{c} \eta\_1 \\ \eta\_2 \end{array} \right) = 0 \tag{38}
$$

with

$$H\_{ii} = \left[\gamma\_i k\_i^2 + \alpha\_i(T)\right] \left[\frac{d}{2} + \frac{\sin(k\_i d)}{2k\_i}\right] \tag{39}$$

and

$$H\_{12} = H\_{21} = -\epsilon\_{12} \left[ \frac{\sin(\frac{k\_1 d}{2} + \frac{k\_2 d}{2})}{k\_1 + k\_2} + \frac{\sin(\frac{k\_1 d}{2} - \frac{k\_2 d}{2})}{k\_1 - k\_2} \right]. \tag{40}$$

The critical temperature of the two-band superconducting film will be determined by the condition

$$H\_{11}H\_{22} - H\_{12}H\_{21} = 0\tag{41}$$

at *T* = *Tc*, which can be explicitly written as

$$\prod\_{i=1,2} \left[ \gamma\_i k\_i^2 + a\_i(T\_c) \right] \left[ \frac{d}{2} + \frac{\sin(k\_i d)}{2k\_i} \right] = \epsilon\_{12}^2 \left[ \frac{\sin(\frac{k\_1 d}{2} + \frac{k\_2 d}{2})}{k\_1 + k\_2} + \frac{\sin(\frac{k\_1 d}{2} - \frac{k\_2 d}{2})}{k\_1 - k\_2} \right]^2. \tag{42}$$

For the two-band superconductor FeSe, we have *Tc*<sup>0</sup> ≈ 8K[5] and the average lattice constant *a* ≈ 0.55 nm [3]. The density of states at the Fermi level for each band are *N*<sup>1</sup> = 0.30 and *N*<sup>2</sup> = 0.65 eV−1, respectively, [8,42]. From the numerical work in reference [43], we can get *λ*<sup>11</sup> = *g*11*N*<sup>1</sup> = 0.30 and the ratio *g*11:*g*12:*g*<sup>22</sup> = 1:0.30:0.37. Then, we have *λ* = 0.054, *λ*max = 0.41 and <sup>12</sup> = 1.1 eV<sup>−</sup>1. With the average Fermi velocities *vF*<sup>1</sup> = 3.7 and *vF*<sup>2</sup> = 4.4 in units of 10<sup>13</sup> nm · <sup>s</sup>−<sup>1</sup> [44], we get *<sup>γ</sup>*<sup>1</sup> = <sup>20</sup> nm<sup>2</sup> · eV−<sup>1</sup> and *<sup>γ</sup>*<sup>2</sup> = <sup>61</sup> nm2 · eV−<sup>1</sup> from Equation (14). Then, from Equation (28), we can obtain the characteristic length scales *A*<sup>11</sup> = (50 nm)−<sup>1</sup> and *A*<sup>22</sup> = (70 nm)−1. For a given film thickness *d*, we first get *ki* from Equation (36). By plugging it into Equation (42), the critical temperature *Tc* as a function of *d* can be calculated numerically and then plotted in Figure 1. It is shown that the critical temperature keeps gradually decreasing with the decreases in *d* and *Tc* ≈ 3 K when the film thickness is reduced to *d* = 3 nm. From Figure 1, we can see that our theoretical results fit the experimental data on Si/SiO2 substrate well.

**Figure 1.** The critical temperature as a function of FeSe film thickness. The experimental data are taken from Ref. [18] (circles) and Ref. [19] (squares) respectively.

At this point, we can also discuss the single-band superconducting system of niobium [45,46], which shows similar behavior to FeSe. The critical temperature of Nb films gradually decreases with a reduction in the thickness from 300 to about 50 nm. A much stronger dependence of the critical temperature on thickness is observed for films thinner than 50 nm. We would also like to apply our theoretical scheme with the similar boundary condition <sup>∇</sup><sup>Ψ</sup> · *<sup>s</sup>*|*<sup>S</sup>* <sup>=</sup> <sup>−</sup>*A*<sup>Ψ</sup> to understand this physical property for singleband superconductor Nb. With the lattice constant *a* = 0.33 nm, the Fermi velocity *vF* = 3.9 in units of 10<sup>13</sup> nm · <sup>s</sup><sup>−</sup>1, the density of states at the Fermi level *<sup>N</sup>* = 1.5 eV−<sup>1</sup> and the dimensionless interaction parameter *λ* = 0.32 for niobium [47,48], we can obtain *<sup>γ</sup>* <sup>=</sup> <sup>7</sup>*ζ*(3)*h*¯ <sup>2</sup>*Nv*<sup>2</sup> *F* <sup>16</sup>*π*2(*kBTc*0)<sup>2</sup> <sup>=</sup> <sup>80</sup> nm<sup>2</sup> · eV−<sup>1</sup> with the bulk critical temperature *Tc*<sup>0</sup> <sup>=</sup> 9.4 K and *A* = *Na γλ* = (52 nm)−1. Thus, with this characteristic length scale, our theoretical scenario can qualitatively explain the strong *Tc* suppression around 50 nm in the single-band superconducting Nb films.

## **4. Conclusions**

In conclusion, we introduced the appropriate boundary conditions in two-band GL theory at the interface between two-gap superconductor and insulator (or vacuum). We also gave a microscopic derivation of these boundary terms based on two-band Bogoliubov–de Gennes formalism. For the two-band superconductor FeSe, we obtained the characteristic length scales of the boundary effect as 50 and 70 nm. The theory can perfectly explain the dramatic suppression of *Tc* when the film thickness is reduced to the same order of these length scales. Our investigation thus suggests that the boundary effect induced by these new terms may play an important role in the research of some iron-based superconducting films.

**Author Contributions:** Conceptualization, C.Y., J.C. and H.H.; methodology, C.Y., J.C. and H.H.; software, C.Y.; validation, C.Y., J.C. and H.H.; formal analysis, C.Y., J.C. and H.H.; investigation, C.Y., J.C. and H.H.; resources, H.H.; data curation, C.Y., J.C. and H.H.; writing—original draft preparation, C.Y.; writing—review and editing, C.Y., J.C. and H.H.; visualization, C.Y.; supervision, J.C. and H.H.; project administration, H.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
