## 1 Introduction  
An underground tunnel is a shallow-buried structure located under the surface of an infinite half-space. Under the action of seismic excitation, it will generate dynamic interaction with the surrounding soil (known as soil-structure interaction, SSI). The methods for solving SSI problems can be roughly divided into two categories: analytical methods and numerical methods. The analytical method mainly refers to the wave function expansion method, which mainly studies regular shape models embedded in a uniform half-space. Although the analytical method is only suitable for relatively simple and regular models, it has an advantage over the numerical method in revealing the essence of the problem. Therefore, enriching the analytical solution library is the goal that scholars have been pursuing. [PERSON] and [PERSON] [1] pioneered the application of the wave function expansion method of electromagnetic waves to the problem of seismic wave scattering and presented the out-of-plane dynamic analytical solution of a single buried tunnel in the whole space. Then, scholars further solved the analytical solutions of the dynamic response of a single tunnel embedded in an isotropic half-space, respectively excited by plane SH-wave, P-wave, and SV-wave [2; 3; 4; 5; 6]. Numerical methods include finite difference method, finite element method, and boundary element method. The numerical method is suitable for arbitrary-shaped tunnels and various site conditions and is more suitable for practical engineering problems. [PERSON] et al. [7] showed the tunnel's response by considering the nonlinear soil conditions. [PERSON] et al. [8] studied the dynamic displacement and stress concentration of an arched tunnel. They found that the dynamic amplifications of the displacements and stresses induced in the tunnel depend crucially on the properties of the surrounding soil, the depth of embedment, and the frequency of the incident simple harmonic disturbance. [PERSON] et al. [9; 10] investigated the three-dimensional dynamic response of an embedded tunnel in a horizontally layered half-space under incident body waves. They extended the comparisons with previous two- and three-dimensional results for the case of a shell embedded in a uniform half-space, and some new numerical results for a shell embedded in a multilayered half-space are presented in a companion paper. [PERSON] et al. [11] presented an analytic solution for the model of a circular lined underground tunnel in elastic half-space excited by in-plane waves, based on the method of complex functions introduced in [12]. The key idea to treat the boundary condition on the ground surface is to find a suitable conformal mapping that transforms the semi-infinite domain of the elastic half-space into a finite circular region, so the original flat half-space surface becomes a closed circle in the image plane which enables the parametrization using azimuthal variable. The advantage of conformal transformation is that some more complicated situations can be considered. Conformal transformation can transform irregularly shaped cross-sections into, for example, circular cross-sections through mathematical transformation, so that the original problem can be solved. The disadvantage of conformal transformation is that the mathematical derivation becomes more complicated. Recently, [PERSON] et al. [13] studied the influence of local topography on tunnel dynamic response. However, all the above studies on soil-tunnel interaction are focused on a single tunnel. Thus, this paper aims to investigate the dynamic interaction between twin flexible tunnels by using the wave function expansion method. The work in this article is an extension of [PERSON]'s research on the seismic response analysis of a single rigid tunnel model [6]. The problem studied in this paper involves the seismic interaction between multiple underground structures, which is currently termed as structure-soil-structure interaction (SSSI) [14; 15; 16; 17; 18; 19].