\section*{Introduction}

While the fifth-generation (5G) wireless network is under deployment globally, researchers have been enthusiastically fostering the future  sixth-generation (6G) wireless network that targets at supporting new and promising Internet-of-everything (IoE) applications ranging from extended reality, automated system, to tactile Internet. These applications impose more stringent performance requirements than 5G, such as ultra-high data rate, global coverage and connectivity, and extremely high reliability and low latency, which may not be fully achieved by existing  technologies for 5G.

To meet the future demands of 6G, an innovative concept of \emph{smart radio environment} has been recently proposed  \cite{Renzo2019Smart,qingqing2019towards}. Specifically, the radio prorogation environment, traditionally deemed to be random and largely uncontrollable, can be dynamically reconfigured  to enhance the wireless communication performance,  by leveraging the digitally-controlled passive metasurface \cite{Renzo2019Smart,qingqing2019towards}.  This promising application of metasurface to wireless communications has rapidly spurred intensive research into the new technology of   \emph{intelligent reflecting surface} (IRS) \cite{qingqing2019towards} or its various equivalents such as  reconfigurable intelligent surface (RIS) \cite{basar2019wireless}.
  The current research on IRS/RIS can be roughly classified into two main categories, namely, hardware design and system modeling \cite{di2020smart}, and communication design and performance study \cite{wu2020intelligent}. {\color{black}Specifically,
IRS can be  utilized to realize a variety of key functions in wireless communication systems by dynamically controlling the amplitudes and/or phases of its reflected signals via its smart controller},  such as  creating (virtual) line-of-sight (LoS) links, improving the channel rank condition,
enhancing desired signal power and/or suppressing co-channel interference by passive beamforming, reshaping the  channel realization and/or statistical distribution, and so on. {\color{black}Moreover, IRS operates in full-duplex (FD) mode with passive reflection only and thus is free of amplification/processing noise as well as self-interference.}
These appealing functions of IRS have motivated substantial studies on applying it to significantly boost the wireless system performance in various scenarios, such as multi-antenna/multi-carrier communications, multi-user
 non-orthogonal multiple access (NOMA), physical-layer security, wireless power transfer (WPT), mobile edge computing (MEC), etc \cite{wu2020intelligent}.

In the existing works on IRS, it is usually deployed at the side of distributed users in the network to help enhance the communication performance with their  aided base stations (BSs),
 as illustrated in the upper part of  Fig.~\ref{Fig:deployment}. In practice, this {\it user-side IRS deployment} is preferred at designated locations/places such as
   hotspot, cell edge, and moving vehicle for enhancing the local coverage. It is worth noting that this deployment strategy is fundamentally due to the   passive signal reflection of IRS  for minimizing its severe

   product-distance path-loss over the two links with its associated remote BS and local user, respectively.  Thus, this deployment is  drastically  different from that for the conventional
    active relay which can amplify  the signal from the source before forwarding it to the destination, and hence

     is usually placed in the middle of them to achieve the best performance.
 Alternatively, the other IRS deployment  strategy that is also able  to minimize  the product-distance path-loss is
  the \emph{BS-side IRS deployment}, shown in the middle part of  Fig.~\ref{Fig:deployment}.
  Note that this IRS deployment resembles that for the conventional reflect-array
  in e.g., satellite communications; while they generally differ in distance with the source (near- versus far-field) as well as main function (low-cost beamforming versus smart channel reconfiguration).
In the following, we compare the user- versus  BS-side IRS deployment strategies in terms of their key metrics including network coverage, reconfigured channel condition, passive beamforming performance, and signaling overhead.
\begin{itemize}
\item[1)] {\bf Network coverage:}
As user-side IRSs serve users in their vicinities only, they are unable to assist users at other locations in the network, especially when the user distributions are unpredictable and vary dynamically. In contrast, the BS-side IRS has the potential to reconfigure channels of all users in the network and thus generally offers  a larger network coverage than user-side IRSs. It is worth noting that for IRSs coated on objects such as walls and facades of buildings, they  can serve users in \emph{half of the space} only (see Fig.~\ref{Fig:deployment}), for both user-/BS-side IRS  deployment strategies.

\item[2)] {\bf Reconfigured channel condition:}
One key function of IRS is  to bypass environmental obstacles for creating virtual  LoS  BS-user  links via IRS smart reflection. To achieve this, IRS can be properly deployed at the user (or BS)  side to establish LoS links with its nearby users (or BS) with high probability owing to their short distances. {\color{black}Besides, the user-side IRS has an additional advantage that, with proper placement, it is more likely to establish  a virtual LoS(-dominant)  link with the BS;}
 whereas it is generally difficult  for a single  BS-side IRS to enjoy  LoS links with all users in the network due to their random locations. Generally speaking, using BS-side or user-side IRSs alone may not be able to guarantee achieving obstacle-free  links between the BS and all users in the network, as shown in Fig.~\ref{Fig:deployment}.
\item[3)] {\bf Passive beamforming performance:}
The BS-side IRS  in general outperforms its user-side counterparts given the same budget on the  total number of reflecting elements  in terms of the maximum passive beamforming gain achievable for each  user. {\color{black}This is because each user can be served by
all reflecting elements at the BS-side IRS,  while it obtains a  smaller passive beamforming gain from its nearby user-side  IRSs only  since the signals reflected by other far-apart user-side IRSs are much weaker due to the higher path-loss \cite{qingqing2019towards}.}
Moreover, the passive beamforming at the BS-side IRS can be more flexibly  adjusted according to all users'  channels and different quality-of-service (QoS) requirements as compared to user-side IRSs that can only cater to their locally served users.
\item[4)] {\color{black}{\bf IRS-BS signaling overhead:}
The IRS reflection design in general requires information exchange between the IRS controller  and its associated BS to tune IRS reflection coefficients. As compared to user-side IRSs, the BS-side IRS generally requires lower signaling overhead with the associated BS, thanks to its much shorter distance  with the BS and thus higher signaling rate, as well as dispensing with the need of coordination among multiple distributed IRSs.}
\end{itemize}

\begin{figure}[t]
\begin{center}
\includegraphics[width=8.7cm]{./hybrid_IRS.pdf}
\caption{Illustration of different IRS deployment strategies.
}
\label{Fig:deployment}
\end{center}
\end{figure}

 \begin{table*}[t!]
\centering
{
\caption{Comparison of three IRS deployment strategies.}\label{Table:Comp}
\resizebox{1\textwidth}{!}{
\begin{tabular}{| c | c |c|c|c|c|}
\hline
\multirow{2}{*}{\textbf{IRS deployment strategy}}&\multirow{1}{*}{\textbf{Network}}&\multirow{1}{*}{\textbf{Reconfigured channel }}&\multirow{1}{*}{\textbf{Passive beamforming}}&\multirow{1}{*}{\textbf{IRS-BS signaling}}&\multirow{1}{*}{\textbf{Double-reflection}}\\
& \textbf{coverage}& \textbf{LoS availability}& \textbf{gain per user}& \textbf{overhead} &\textbf{link }\\
\hline
User-side deployment & Relatively small & High& Relatively small&High & No \\
\hline
BS-side deployment & Large &  Relatively low& Large & Low& No \\
\hline
Hybrid deployment & Very large & Very high& Very large & High& Yes\\
\hline
\end{tabular}}
}
\end{table*}

Motived by the above, we further propose a more general
 \emph{hybrid IRS deployment} strategy by combining the complementary advantages of the BS- and user-side IRSs with more flexibility to trade-off between them, as shown in the lower part of Fig.~\ref{Fig:deployment}.

Moreover,  besides the  IRS single-reflection (i.e., BS-IRS-user) links in the user-/BS-side deployment cases,
the proposed hybrid IRS deployment leads to the additional inter-IRS reflection link between the BS-side IRS and each of the user-side IRSs (e.g., see the  BS-IRS $1$-IRS ${2}$-user link shown in Fig.~\ref{Fig:deployment}).
The \emph{double-reflection} links can be exploited to establish more available LoS paths between the BS and its served users, especially when
both the direct and single-reflection links between them  are severely  blocked.
 More interesting, it has been shown in  \cite{Han2020Cooperative} that under the LoS channel setup,
 the double-reflection link can achieve a much higher asymptotic passive beamforming gain as compared to the single-reflection link given the same total number of reflecting elements, denoted by $N$, as $N$ becomes sufficiently large (i.e., in the order of $\mathcal{O}({N^4})$ versus $\mathcal{O}(N^2)$), despite the fact that the former suffers more path-loss than the latter due to the double (versus single) signal reflection.

Although  the hybrid IRS deployment potentially offers better flexibility and hence superior   performance as compared to the user-/BS-side IRS deployment alone, it also incurs higher complexity in design and implementation, explained as follows.
First, how to optimally assign the reflecting elements at the BS-side IRS to serve users directly (i.e., via the BS-side single-reflection links only) or with some of the user-side IRSs cooperatively (i.e., via the double-reflection links over them) is an important but difficult problem to solve. This problem is further complicated with the joint design of the  active beamforming of the BS and the passive beamforming of all IRSs.
Second,
 how to efficiently acquire the necessary channel state information (CSI) for the above-mentioned IRS association and joint active/passive beamforming  design is more  practically   challenging than the case with user-/BS-side IRSs only,

  due to the presence of both the single- and double-reflection links and their intricate coupling. Third, both the BS- and user-side IRSs need to be properly placed in the network, and the reflecting elements allocation among them should be optimized for maximally enhancing the network performance, which is also a hard problem to solve in practice.  In  Table~\ref{Table:Comp}, we summarize the comparison of the three IRS deployment strategies  in various main  aspects.

\section*{Design Challenges and Promising Solutions }
 In this section, we present the main challenges in designing the hybrid IRS deployment, including the IRS-user  association and mode selection, CSI acquisition and beamforming optimization, as well as  IRS placement and elements allocation. Promising approaches are also proposed to solve them.

\section*{\small{IRS-User  Association and Mode Selection}}

\begin{figure*}[t!]
\begin{center}
\includegraphics[width=14.2cm]{./mode.pdf}
\caption{Illustration of different  communication modes in an IRS-aided system.}
\label{Fig:mode}
\end{center}
\end{figure*}

Under the hybrid IRS deployment, the foremost challenge in optimizing the wireless network performance is the design of  IRS-user  association. Compared to the user-/BS-side IRS deployment, IRS-user  association under the hybrid IRS deployment is  more involved. In this case, each of the users may be associated with no IRS, the BS-side IRS only, some user-side IRS only, or both the BS- and user-side IRSs, depending on its communication mode, which can be direct communication with the BS only, communication with the BS over the user-/BS-side IRS single-reflection link only,  or over both the (two) single-reflection and (one) double-reflection links together, as illustrated  in Fig.~\ref{Fig:mode}.
 Generally speaking, the design of IRS-user  association and (user) mode selection (IA-MS) needs to take into account various practical factors

 such as users' QoS requirements, channel conditions, co-channel interference levels, and locations in the network. Besides, the IA-MS design is also coupled with
 the active beamforming design at the BS and other communication resources (e.g., power, bandwidth) allocation,  users' multiple access scheme and transmission scheduling, as well as the reflection coefficients at all IRSs in the network.

  Moreover, as the users move within the same cell, their associated BS may be unchanged, while their associated IRSs and communication modes usually  need to be adaptively adjusted, especially for high-mobility users without any user-side IRS attached.

However, the optimal IA-MS design is practically challenging and may be even impossible. This is because 1) it requires the full  CSI of all involved communication links in the network which is difficult to obtain in practice; and 2) it is intricately coupled with the active/passive  beamforming design and hence cannot be solved alone, even with perfect CSI. To tackle the above difficulties, a practically appealing approach is to   leverage the \emph{statistical CSI} for designing the
 IA-MS based on a practical  \emph{surface-partition} scheme. Specifically,
for the  BS-side IRS, the statistical CSI (e.g., the deterministic large-scale fading channel components  and fading channel spatial correlation matrix) of its links with the associated BS and other user-side IRSs (as well as their locally served users due to close proximity)    can be first obtained (by e.g., equipping the IRS with low-cost sensors). Moreover, the BS-side IRS can be partitioned into multiple subsurfaces, which consist of equal or unequal number of reflecting elements (see Modes 3 and 4 in  Fig.~\ref{Fig:mode}). As such,
  instead of associating users to  the entire BS-side IRS with/without any helping user-side IRSs at the same time,
a low-complexity scheme is to assign each subsurface to one  user only at each time, based on the statistical CSI, users' QoS requirements, and their estimated link signal-to-noise ratios (SNRs) with/without a helping user-side IRS. {\color{black}{For example, under the binary LoS/non-LoS(NLoS) channel model and based on the large-scale channel gains (e.g., path loss) for all involved links, the received signal power of each user's link with respect to the number of IRS subsurfaces assigned  can be estimated. Based on them for all users, an integer programming can be formulated to optimize the IRS surface partition for optimizing the performance, e.g.,  maximizing the minimum received signal power among all users. In this case, the user with a smaller channel power gain should be allocated with more IRS subsurfaces.}}
While the surface-partition scheme also applies to each of the user-side IRSs, its surface size is practically much smaller than that of the BS-side IRS, thus rendering its channels with locally served users more likely to be spatially  correlated.
 Therefore, instead of applying the surface-partition scheme,  an efficient IRS-user  association method for the user-side IRS is to leverage its time-varying reflections to serve its locally associated users {\color{black}{over different time.}}
 Accordingly, the SNRs of these users can be estimated and then fed back to the BS for facilitating the design of the BS-side IRS association as well as determining whether there is a need for each user to be associated with its nearest user-side IRS, if any.

It is worth noting that although the proposed IA-MS design based on statistical CSI and surface partition  is suboptimal in general, it significantly reduces the design complexity and also helps  simplify the subsequent CSI acquisition and  beamforming design  (as will be detailed in the next subsection).
 Moreover, the proposed IA-MS design can be conveniently adjusted over time for different users  according to their movement and/or sporadic traffic by simply re-assigning each IRS subsurface to a different user, without affecting other subsurfaces and their user associations.

Despite the above advantages, the proposed  IA-MS design also faces several new challenges. First, given user locations, how to solve the IRS-user association problem based on statistical CSI is practically challenging, due to the co-existence of both the BS- and user-side IRSs that introduce complex channels (with both single- and double-{\color{black}reflection} links), as well as multiple user communication modes.
 Besides, the performance loss of the suboptimal~surface-partition approach needs to be characterized as compared to the optimal IRS-user associations without surface partition.

 Moreover, the IA-MS design based on statistical CSI can be further improved/refined based on the actual SNRs of the users with the active/passive beamforming and other communication resources allocation implemented. However, how to devise an efficient on-line algorithm to achieve this goal is still open.

\section*{\small{CSI Acquisition and Beamforming Optimization}}

With the proposed  IA-MS design, all users' IRS associations and communication modes are determined. Then, efficient CSI acquisition and joint active/passive  beamforming designs  can be implemented  for optimizing the users' QoS performance,  which can be roughly classified into two categories, namely, \emph{channel-estimation based} and \emph{beam-training based} methods.

The channel-estimation based method needs to firstly estimate the CSI of all involved  links and then design the joint active/passive  beamforming based on estimated CSI. For IRS channel estimation under the user-side IRS deployment, there are two main approaches in the existing literature based on different IRS configurations, depending on whether sensing devices  are integrated into IRS or not. First, for semi-passive IRS mounted with sensing devices for receiving signals, {\color{black}the CSI of  the BS/user$\to$IRS links  can be constructed from that of the BS/user$\to$sensor links by using e.g., data-interpolation and  channel-calibration techniques.} However, this method only works for time-division duplexing (TDD) systems based on channel reciprocity, while it is infeasible for frequency-division duplexing (FDD) systems since the CSI of the  links from IRS to the BS/users  cannot be obtained with sensing devices/receivers  only. Second, for (fully) passive IRS without sensing devices integrated, a viable channel estimation approach is by estimating the cascaded BS-IRS-user reflection channel at the BS/users based on uplink/downlink pilot signals, respectively, while the IRS should dynamically tune its training reflection over time to facilitate the cascaded channel estimation \cite{wu2020intelligent}. However, when extending this method to the general case with the hybrid IRS deployment, a new  challenge
arises due to  the double-reflection link since its cascaded CSI is more difficult to estimate than that of the user-/BS-side IRS single-reflection links.  This is because not only more channel coefficients are involved in the double-reflection link, but also the double-reflection link is intricately coupled with the two single-reflection links (see Mode 4 in Fig.~\ref{Fig:mode}). To address this problem, efficient channel estimation schemes have been recently proposed in \cite{you2020wireless} and \cite{zheng2020efficient} to estimate the double-reflection channel for the single-user case, and both single- and double-reflection channels for a cluster of users, respectively.
  However, how to generalize the results in the above works to the case with arbitrary user  channels and locations  under the hybrid IRS deployment  still needs further investigation. {\color{black}Moreover, it is also interesting to study how to apply deep learning techniques to improve the channel estimation performance and training efficiency under the hybrid IRS deployment \cite{liu2020deep}.}
 Based on the estimated CSI, the active beamforming at the BS and passive beamforming at all IRSs can be jointly designed \cite{Zheng2020DoubleIRS}. {\color{black}Specifically,  it has been shown that deploying two cooperative IRSs near the BS and users achieves superior performance to the case with single-IRS deployment in terms of the maximum SNR  and multi-user effective channel rank \cite{Zheng2020DoubleIRS}.    However, how to devise

robust beamforming against CSI estimation errors  under the general hybrid IRS deployment is still an unsolved problem in the literature \cite{zhou2020framework}.}

Alternatively, the beam-training based method does not require the CSI explicitly. Instead, it is a codebook-based beamforming method that directly searches for the best beam in the spatial domain that yields the optimum performance, which is thus of lower complexity and more scalable with the increasing number of IRS reflecting elements in the network as compared to the channel-estimation based methods {\color{black}{due to the finite number of beam codewords}}. {\color{black}While an efficient IRS beam training scheme has been recently proposed in \cite{wang2021joint} that is  applicable to systems under the user-/BS-side  IRS deployment, its extension to the hybrid IRS deployment is still open, which is more challenging due to the additional inter-IRS channels and thus calls for new beam training designs. For example, to reduce the beam training overhead, an efficient approach is by designing the hybrid offline-online beam training, where the beam training for quasi-static BS-IRS and inter-IRS channels is performed offline prior to data transmissions, while the IRS beam directions to the mobile users are trained in real time.}
{\color{black} In practice, the beam training method can be used to generate a proper fingerprinting database that collects a set of optimal beamforming designs at some user positions, with which
deep learning techniques can be leveraged to predict the mapping between the position of an arbitrary user and its corresponding optimal beamforming.} More research along this direction is needed for further exploration.

\section*{\small{IRS Placement and Elements Allocation}}

While IRS deployment is assumed to be  fixed in the preceding subsections, we consider in this subsection the  design of IRS deployment, including the IRSs' placement and their reflecting elements' allocation.

First, for each user-/BS-side  IRS,  there are generally two placement strategies, namely, the \emph{centralized placement} where all its reflecting elements are placed on one single surface, versus the \emph{distributed placement} where the available reflecting elements  form multiple smaller-size sub-IRSs that are placed  at different locations around the BS/users. Compared with centralized IRS placement, distributed IRS placement in general achieves better
 channel conditions, e.g., increasing the LoS probabilities of the IRS-user and/or the BS-IRS  links, and improving their  channel ranks by exploiting the spatial path diversity.
   However, distributed IRS placement generally requires higher signaling overhead and hardware cost due to more IRS controllers that need to be  installed, each for one sub-IRS. A thorough performance comparison of the above two placement strategies under the general hybrid IRS deployment is still lacking in the literature, which thus calls for future studies.

Next, under  the hybrid IRS deployment with centralized/distributed IRS placement, a unique design challenge lies in how to optimally allocate the  available reflecting elements  between the BS- and user-side IRSs for optimizing the network performance {\color{black}given a fixed number of total reflecting elements}.  {\color{black}To draw useful insights into this elements-allocation problem, we first consider two simplified system setups.} First, for a multi-user system where  all reflecting elements are assumed to be deployed  near either  the BS or distributed users, it has been shown in \cite{zhang2020intelligentcapa} that the BS-side IRS deployment yields superior rate performance to the user-side deployment under a ``twin" channel setup (i.e., symmetric  channel realizations for fair comparison). Second, for a single-user system,
 it was revealed in \cite{Zheng2020DoubleIRS} that, the available reflecting elements should be split to form two comparable-size IRSs, each deployed near the BS and the user, respectively, for reaping the double-reflection multiplicative beamforming gain.

The above  two results indicate that for rate maximization under the hybrid IRS deployment, the available reflecting elements should be properly allocated to both the BS- and user-side IRSs.  {\color{black}However, under the proposed hybrid IRS deployment with generally multiple users and various channel conditions (e.g., the LoS availability) in the links between the users/BS and BS-/user-side IRS,
how to solve this general IRS elements allocation problem is still open. A viable approach to solve this hard problem is by firstly determining the IRS-user association based on the statistical CSI and then optimizing the IRS elements allocation by using the  continuous variable  relaxation and integer-rounding techniques for e.g., maximizing the minimum rate among all users.
  Moreover, for a large-scale network with massive randomly distributed users, the information-theoretical approach adopted in \cite{Zheng2020DoubleIRS} and \cite{zhang2020intelligentcapa}  for optimizing IRS elements  allocation might be computationally formidable. As such, new mathematical tools such as stochastic geometry and deep learning may need to be invoked to facilitate the network-level analysis and efficient elements allocation.} The study on IRS placement and elements allocation is still at an initial stage and thus deserves more research efforts  in the future.

\section*{Numerical Results}

Numerical results are provided in this section to demonstrate the effectiveness of the proposed hybrid IRS deployment.
We consider a typical IRS-aided single-cell system under the general hybrid IRS deployment as shown in Fig.~\ref{system_hybridIRS}, where one BS-side IRS (denoted as IRS $1$) of $N_1$ reflecting elements  and two user-side IRSs (denoted as IRS $2$ and $3$) of $N_2$ and $N_3$ reflecting elements, respectively, are deployed to assist in the uplink transmission from three single-antenna users to  a $40$-antenna BS,  {\color{black} with the locations of all nodes specified in Fig.~\ref{system_hybridIRS}}. Note that the considered hybrid IRS deployment reduces to the user-side and BS-side deployment,  respectively,  when $N_1=0$ and $N_2=N_3=0$. As for  IRS coverage, the BS-side IRS $1$ can serve the cell-edge users $1$ and $2$ in its serving half-space, while user $3$ is outside its serving coverage. On the other hand,  user-side IRSs $2$ and $3$ can serve users $2$ and $3$, respectively, in their local coverage.
{\color{black}Moreover,  we assume that the links between the BS/users and their nearby serving IRSs follow the LoS channel model with path-loss exponent of $2.2$; the links between the BS/users and their far-apart IRSs are modeled by the LoS-dominant Rician fading channel with Rician factor of $5$~dB and path-loss exponent of $2.5$; and the BS-user links follow the Rayleigh fading channel model with path-loss exponent of $3$.
All the users are assumed to have equal transmit power of $5$~dBm and the noise power at the BS receiver is  $-75$~dBm.}

\begin{figure}[!t]
	\centering
	\includegraphics[width=7cm]{system_hybridIRS_v4.pdf}

	\caption{System setup for the simulation (top view).}
	\label{system_hybridIRS}
\end{figure}

\begin{figure}[!t]
	\centering
	\includegraphics[width=9cm]{hybridIRS_compV2.pdf}

	\caption{Achievable rates of the three typical users as well as their min-rate under different IRS deployment strategies.}
	\label{hybridIRS_compV2}
\end{figure}

Let $N=400$ denote the budget on the total number of reflecting elements for all IRSs, i.e., $N_1+N_2+N_3=400$. We consider the following elements allocation for the three IRS deployment strategies under the considered system setup: 1)   $N_2=N_3=N/2=200$ for the user-side IRS deployment (only); 2) $N_1=400$ for the BS-side IRS deployment (only); and 3) $N_1=300$ and $N_2=N_3=50$ for the hybrid IRS deployment. Moreover, for the BS-side IRS, we apply the exhaustive search to find its optimal surface partition {\color{black}{for maximizing the minimum achievable rate (min-rate) among all users}}. In addition, the multi-antenna BS applies the maximal-ratio combining (MRC) receive beamforming  to the received signal from  each user.

In Fig.~\ref{hybridIRS_compV2}, we compare the min-rate  in bits per second per Hertz (bps/Hz) among the three users for achieving user fairness  under the three IRS deployment strategies.
First, it is observed that under the user-side IRS deployment, user $1$ achieves a significantly lower  rate than users $2$ and $3$, thus becoming the min-rate performance bottleneck in the system. This is expected since the cell-edge user $1$ is outside the local coverage of any user-side IRS and thus served by the BS directly. Next,  the BS-side IRS deployment strategy is observed to effectively enhance the achievable rate of user~$1$ without sacrificing much the rate performance of user $2$, thanks to its great flexibility in the surface partition and passive beamforming.  In this case, user $3$ becomes the min-rate performance bottleneck since it is not in the serving half-space of the  BS-side IRS (i.e., IRS $1$) and thus cannot enjoy the passive beamforming gain. Third, it is observed that with proper  reflecting elements allocation as well as efficient BS-side IRS surface partition, the hybrid IRS deployment strategy achieves  a much higher min-rate than both the user- and BS-side IRS deployment strategies due to the following reasons. First, in the considered setup, all the users can be covered by at least one (BS-side or user-side) IRS under the hybrid IRS deployment, which demonstrates its more  effectiveness in terms of  network coverage. Second, user $2$ has the potential to enjoy the higher-order passive beamforming gain from the double-reflection link (as compared to the single-reflection link), and hence more reflecting elements at the BS-side IRS can be assigned to associate with user $1$ for improving its achievable rate over the  single-reflection link only, thus more  flexibly balancing their rates for improving the system min-rate.

\section*{Conclusions}
In this article, we provide an overview of the IRS deployment design in wireless networks from a communication perspective.
We first
compare  the conventional   BS- versus  user-side IRS deployment strategies, and then introduce a new   hybrid IRS deployment strategy by combining their complementary advantages.
The main challenges for designing the IRS-aided wireless network under the general hybrid deployment are discussed, including the IRS-user association and mode selection, IRS channel acquisition and beamforming optimization, and IRS placement and elements allocation. Promising methods to tackle these challenges as well as directions for future investigation are provided.
 Numerical results are also presented  to corroborate our discussions and show the performance gains of the hybrid IRS deployment as compared to its one-sided counterparts.   In future work,
it is also interesting to study
the more general multi-IRS aided communication system with multi-hop (i.e., more than two hops) signal reflection. In this case, IRS beamforming needs to be jointly designed with data routing over multiple selected IRSs for serving each user  \cite{mei2020cooperative},  as illustrated in Fig.~\ref{Fig:routing}. More research is  needed to develop advanced techniques for designing the multi-IRS cooperative multi-beam multi-hop routing for multiple users with joint active/passive beamforming optimization, so as  to realize the full vision of smart radio environment in future wireless systems. {\color{black}Moreover, in this article, we mainly focus on the IRS deployment design from the signal processing and optimization perspective, while other useful tools such as stochastic geometry and machine learning \cite{liu2020ris} can also be applied to obtain efficient IRS deployment designs, especially for large-scale wireless networks, which deserve further investigation.}

\begin{figure}[t!]
\begin{center}
\includegraphics[width=8.5cm]{./beam_routing.pdf}
\caption{Multi-beam multi-hop routing in the multi-IRS aided multi-user communication system via joint active/passive beamforming.}
\label{Fig:routing}
\end{center}
\vspace{-10pt}
\end{figure}


\section{Introduction}

Interactions between shock waves and turbulent flows appear in multiple high-speed applications including supersonic flight and combustion processes \citep{Andreopoulos2000}.
Specifically, in shock wave boundary layer interactions (SBLI), the imposed adverse pressure gradient over the wall induces flow separation, which is associated with a reflected shock presenting highly unsteady low-frequency fluctuations.
This situation occurs in many practical configurations such as turbomachinery flows, transonic buffeting, overexpanded rocket nozzles and supersonic inlet isolators.
Besides affecting the flow downstream of the interaction zone, these oscillations are responsible for large pressure levels and severe structural loading.
The source of these broadband shock fluctuations has been the subject of extensive debate and comprehensive reviews over the past decades \citep{Delery1985,Dolling2001,Dussauge2006,Babinsky2011,Clemens2014,Gaitonde2015}.

Two mechanisms are recognized to induce reflected shock unsteadiness: upstream boundary-layer forcing and downstream feedback \citep{Clemens2014}.
Among numerous investigations \citep{Plotkin1975, Brusniak1994, Beresh2002}, a strong correlation between the upstream boundary layer streamwise velocity and shock motion was found in the experiments of \citet{Ganapathisubramani2009}.
Three-dimensional measurements from \citet{Humble2009} demonstrated that low and high-momentum, streamwise-elongated zones convected within the incoming boundary layer affect directly the reflected shock structure.
On the other hand, there is compelling evidence for the influence of the downstream separated flow on the shock motion.
The large-eddy simulations (LES) and linear stability analysis of \citet{Touber2009} suggested the existence of unstable modes and upstream traveling waves associated with a global instability mechanism to be responsible for the low-frequency shock oscillations.
These oscillations were found to be highly correlated to the reattachment pressure measurements performed by \citet{Dupont2006}.
The direct numerical simulation (DNS) of \citet{Wu2008} also pointed to a strong correlation where the shock motion lags the reattachment fluctuations.
The mass-entrainment model proposed by \citet{Piponniau2009} scales the low-frequency content with the separated flow region and highlitghts the importance of the downstream flow on selecting the shock unsteadiness. It is worth noting that a low-frequency feedback from the reattachment region up to the vicinity of the reflected shock was documented for transitional SBLI \cite{sansica_sandham_hu_2016,larchevequeaiaa}. Although the flow features listed above could suggest that a similar mechanism may also be at play for turbulent SBLI, this seems not to have been demonstrated in the current literature.

The low-frequency unsteadiness and the associated mechanisms discussed above were observed for impinging oblique shocks and compression ramp configurations, in both numerical and experimental investigations.
Equally important is the spatial organization of the low-frequency motion in the vicinity of the reflected shock and separation bubble.
The low-pass filtered velocity fields from the simulations of \citet{Priebe2012} highlighted a strong quasi two-dimensional behavior of the reflected shock motion using spanwise-averaging.
The dynamic mode decomposition analysis (DMD) conducted by \citet{Priebe2016} shows that the shock unsteadiness was reproduced using a number of low-frequency DMD modes.
Their spatial imprint qualitatively agreed with the linear unstable modes found in \citet{Touber2009}.
These modes exhibit streamwise-elongated structures in the downstream separation bubble of the compression ramp.
Similar streamwise Görtler-like vortices were found in the DMD analysis from \citet{Pasquariello2017} for a strong impinging shock wave.
The global stability analysis and DMD computed over spanwise-averaged velocity fields from \citet{Nichols2016} shows a striking resemblance with the recirculation bubble breathing during one cycle of the oscillation.

In the present work, we exploit the LES database from \citet{Jiang2017} where an impinging shock wave interacts with a turbulent boundary layer at Mach number $M=2$ and $Re_\theta\simeq5000$.
This database spans a time domain an order of magnitude larger than the current existing DNS or LES simulations investigated for this configuration \citep{Pirozzoli2006,Touber2009,Pasquariello2017}, corresponding roughly to more than 300 low-frequency cycles.
First, we present a causal analysis in the framework of linear transfer functions to identify the link between the low-frequency shock motion, upstream and downstream mechanisms.
To shed further light on the spatial structure of the low-frequency unsteadiness, we perform a spectral proper orthogonal decomposition (SPOD) of the spanwise-averaged velocity fields in the interaction zone.
To the authors' knowledge, little is known about this optimal spatial-spectral data decomposition for SBLI configurations.
In section II, a description of the numerical database is provided accounting for the whole set of simulation and flow parameters.
The key spectral features of the shock motion are analyzed in section III.
The mathematical framework to compute the linear transfer functions, causal analysis and spectral decomposition modes is presented in section IV.
We report the prediction and causality results in section V, and the spectral eigenvalues and modes in section VI.
Finally, the concluding remarks are described in section VII.

\section{Numerical data base}

The current large-eddy simulation (LES) reproduces the shock wave/turbulent boundary layer interaction investigated in the experiments of \citet{Schreyer2015}.
The full description of the numerical simulation was presented in \citet{Jiang2017}, together with a thorough validation against wind-tunnel measurements.
Here, we provide a brief overview of the configuration and main flow parameters.

The numerical setup mimics the experimental shock generated by a flow deflection of $8.5^\circ$ impinging on a $M_{\infty}=2$ turbulent boundary layer with $Re_\theta\simeq 5000$.
The complete set of flow parameters are listed in table~\ref{tab:parameters}.
The incident shock corresponds to a pressure ratio of 1.58 with an angle of approximately $38^\circ$, as indicated in figure~\ref{fig:Mean_Mach}.
All simulations were performed using the FLU3M codes from ONERA.
The code relies on a finite volume discretization in space and an implicit time integration, both being second--order accurate.
The space scheme is designed to minimize the numerical dissipation by adding the dissipative part of the Roe scheme to a centered scheme in regions with strong compressibility/low vorticity, as identified by means of Ducros'sensor \citep{Roe1981,Ducros1999}.
The time integration is performed with a maximum CFL number of 11, using 7 sub-iterations to solve the non-linear system in order to ensure that residuals are reduced by at least 1.5 order of magnitude.
The LES modeling is built from an implicit grid filtering coupled with an explicit subgrid modeling through the selective mixed-scale model that was successfully used in previous studies of shock wave/boundary layer interactions \citep{Garnier2002,Agostini2015}.

\begin{figure}
\centering
\includegraphics[scale=0.65]{./Sections/Mean_Mach_Sonic_Line_Normalized_Final.pdf}
\caption{Mean Mach number in the streamwise ($x$) wall-normal ($y$) longitudinal plane at $z=0$. The turbulent boundary layer (TBL) thickness ($\delta_{99}$) is shown in the first region of the simulation, and the sonic line is represented in the interaction zone.}
\label{fig:Mean_Mach}
\end{figure}

\setlength{\tabcolsep}{11pt}

\begin{table}
\begin{center}
\caption{Flow parameters.}
\begin{tabular}{cccccccc}
\\
$M_{\infty}$ & $p_o$ & $T_o$ & $U_{\infty}$ & $\delta_{99}$ & $Re_{\theta}$ & $L$ & $W$   \\
\hline \\
2 & $40500\,\text{Pa}$ & $295\,\text{K}$ & $510\,\text{ms}^{-1}$ & $11\,\text{mm}$ & 4850 & $57\,\text{mm}$ & $33\,\text{mm}$ \\
\hline
\end{tabular}
\label{tab:parameters}
\end{center}
\end{table}

The wall is modeled as adiabatic.
Non-reflective boundary conditions relying on characteristic variation in space are set at the inflow, outflow and upper boundaries, while periodicity is used in the spanwise $z$ direction.
The incident shock is created by enforcing the pre/post-- shock states on the upper boundaries, following \citet{Garnier2002}.
A fully turbulent inflow condition is obtained by adding stochastic velocity fluctuations to the mean profiles by means of a Synthetic Eddy Method~(SEM) \citep{Jarrin2006}.
Temperature and density fluctuations are obtained from the velocity fluctuations assuming strong Reynolds analogy and linearised ideal gas law with negligible pressure fluctuations \citep{Agostini2012}.
The space and time scales required by SEM are obtained from a pre-existing LES of the same $M=2$ boundary layer without shock impingement that encompasses the laminar-to-turbulent transition process \citep{Schreyer2016}.
Comparisons between the two LES demonstrate that a relaxing length of about 10 boundary layer thicknesses $\delta_{99}$ is sufficient to achieve fully turbulent canonical first-- and second--order statistics from the inflow boundary condition.

The reference mesh was designed by setting the inflow location $16\delta_{99}$ upstream of the shock impingement location.
The grid resolution and spanwise extent $W$ of the computational domain were increased with respect to previous studies \citep{Agostini2015}.
The corresponding resolution in wall units were $\Delta x^+=28$~(streamwise), $\Delta y_{min}^+=0.85$~(wall--normal) and $\Delta z^+=12$~(spanwise) for $W=3\delta_{99}$.
Two additional meshes were tested for a grid-convergence study: a refined one with resolution increased by $40\%$, $20\%$ and $30\%$ in the streamwise, wall--normal and spanwise directions, respectively, and an enlarged spanwise domain corresponding to $W=6\delta_{99}$.

A key characteristic of the numerical database is that data is sampled from the reference mesh over approximately 300 periods of the typical low-frequency unsteadiness, as will be discussed later.
It allows the computation of statistical and spectral quantities, and even conditional statistics with a very low level of statistical uncertainty.
Computations from the two additional meshes were performed for 10 to 18 low-frequency cycles.
These tests were long enough to achieve statistical convergence for validation purposes.

A complete description with validation results from the experimental database and grid convergence can be found in \citet{Jiang2017}.
All the computations and experiments have an interaction length of $L = 5.2\delta_{99}$ and the mean velocity profiles upstream of and within the interaction region agree very well.
Downstream of the interaction region, the LES velocity fields relax toward a canonical boundary layer profile slightly more slowly than in the experiments.
This was explained by the finite width of the experimental setup with side walls compared to the periodicity set in the computations.

Second-order velocity statistics also show very good agreement between the three computations and the particle image velocimetry (PIV) measurements, with the exception of the streamwise fluctuations in the incoming boundary layer and the initial stage of the mixing layer developing over the separated region.
In the latter, the LES underestimate the $\overline{u'u'}$ level by 40\% whatever the mesh used.
Comparisons between the streamwise velocity spectra obtained from the LES and the experiments show that the underestimation originates solely in the low--frequency band $f\,\delta_{99}/U_\infty<0.15$.
It is therefore postulated that it was due to the intrinsic inability of the present SEM method to generate large-scale convective structures found in the experiments, possibly originating in the nozzle upstream of the flat surface over which the boundary layer develops.
The otherwise very good agreement between LES and experiments up to the spectral quantities in multiple locations indicate that for the present flow parameters such fluctuations have very limited impact on the dynamics of the interaction.

\section{Spectral analysis of the shock motion}

A space-time spectral analysis is considered to identify the key features of the shock motion.
The variables are normalized by the interaction length scale $L$ and the free-stream velocity $U_{\infty}$.
The dimensionless time is defined as $t=t^* U_{\infty}/L$, while the normalized frequency reads $St=fL/{U_{\infty}}$, corresponding to the Strouhal number.

First, we assess the dominant frequency in the vicinity of the shock.
The streamwise evolution of the pre-multiplied spectrum of pressure fluctuations is shown in figure \ref{Spectra}.
The spectra is defined as $S_{pp}=\langle \hat{p}\hat{p}^{*}\rangle$, where $\hat{.}$ indicates the Fourier transform from time to frequency $f$, and $\langle{.}\rangle$ defines an ensemble averaging.
In the upstream boundary layer ($x<-0.25$), the spectrum presents mainly high-frequency content ($St>1$) linked to the incoming turbulent eddies.
However, in the vicinity of the shock position ($-0.2<x<0$), a low-frequency broadband range emerges, and is centered approximately at $St=0.03$.
The large gap between these frequency scales was reported in previous investigations \citep{Touber2011,Dupont2006}.
As discussed later in the manuscript, a low-pass filter will be applied to the data to isolate these shock oscillations.

\begin{figure}
\begin{center}
\includegraphics[scale=0.2]{./Sections/Spectra.pdf}
\end{center}
   \caption{Streamwise evolution of the pre-multiplied and normalized spectrum of the wall pressure fluctuations $p$ ($z=0$).}
  \label{Spectra}
\end{figure}

To shed further light on the spatial structure of these oscillations, the space-time spectrum is computed for the streamwise velocity fluctuations.
Here, $S_{uu}(St,\beta)=\langle \hat{\hat{u}}\hat{\hat{u}}^{*}\rangle$, where the double hat indicates a double Fourier transform, from time to $f$ and from $z$ to the spanwise normalized wavenumber $\beta$.
Figure \ref{2DSpectra} presents the two-dimensional spectra in the near-wall region at four streamwise locations $x = -0.35$, 0.0, 0.2 and 1.4.
These locations correspond to the upstream boundary layer, the approximate shock position, the recirculating bubble and downstream reattachement, respectively.
Upstream of the shock, the flow is dominated by broadband turbulent near-wall streaks.
In the approximate shock position, the fluctuations are quasi two-dimensional ($\beta=0$).
As discussed above, this region is predominantly governed by low-frequency fluctuations, which can therefore be isolated via spanwise averaging.
This property will be considered and analyzed in our flow predictions.
The values of $\beta$ increase in the recirculating flow region and further downstream.
The velocity fluctuations exhibit mostly a three-dimensional behavior, and the characteristic frequencies are also higher in this region.
The simultaneous increase of both $St$ and $\beta$ is related to the development of shear-layer vortices and their break-up into 3D eddies.
This behavior can be also attributed to the existence of G\"ortler-type structures developing along the boundary-layer \cite{loginov_adams_zheltovodov_2006,Pirozzoli2006}.

\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.45\textwidth]{./Figs/Updated_Figures/Spectrum_figures/Spectra_2D_u_x035577-eps-converted-to.pdf}}
\subfigure[]{\includegraphics[width=0.45\textwidth]{./Figs/Updated_Figures/Spectrum_figures/Spectra_2D_u_x0002727-eps-converted-to.pdf}}\\
\subfigure[]{\includegraphics[width=0.45\textwidth]{./Figs/Updated_Figures/Spectrum_figures/Spectra_2D_u_x021782-eps-converted-to.pdf}}
\subfigure[]{\includegraphics[width=0.45\textwidth]{./Figs/Updated_Figures/Spectrum_figures/Spectra_2D_u_x14188-eps-converted-to.pdf}} \\
\end{center}
   \caption{Two-dimensional spectra at $x=-0.35$, 0.0, 0.2 and 1.4, respectively in (a) to (d).}
  \label{2DSpectra}
\end{figure}

\section{Methods}

\subsection{Linear transfer function: single-input single-output analysis}

A linear, frequency-domain identification method is applied to relate an input signal $I(t)$ to an output $O(t)$ separated in the streamwise/wall-normal plane.
The supporting assumption for the method is the existence of a linear function that maps the input to the output locations.
A similar approach was applied to estimate velocity fluctuations in turbulent jets and boundary layers \citep{sasaki2017real,sasaki2019transfer}.
Our goal is to identify the causal mechanisms behind the low-frequency shock unsteadiness.
The problem is formulated in the frequency domain,
where the optimal frequency response reads

\begin{equation}
\label{tfequation1}
G_{IO}(f)=\frac{S_{IO}(f)}{S_{II}(f)},
\end{equation}

\noindent where $S_{II}(f)$ and $S_{IO}(f)$ are the auto and cross spectra of the input and output signals, respectively, and $G_{IO}(f)$ refers to the empirical transfer function computed in the frequency domain \citep{bendat2011random}.
These quantities were calculated from the expected values of $\hat{I}(f)\hat{I}^{*}(f)$ and $\hat{I}(f)\hat{O}^{*}(f)$, obtained from an ensemble averaging, via the Welch's method. The spectra was obtained using 230 blocks of 2000 elements each with an overlap of 75\%, which results in a frequency resolution of 0.015. Similar results were obtained using up to 5000 elements in each block.

Equation \eqref{tfequation1} is referred to as the $H_1$ estimator of the system, and it minimizes the error due to measurement noise in the output \cite{rocklin1985comparison}.
Other formulations, such as the $H_2$ or $H_{\nu}$ estimators, exhibit different performances in terms of sensor noise minimization.
However, they are expected to perform equally well for this type of estimation, which does not consider the presence of measurement uncertainties and is based in simulation data.
The $H_1$ estimator presents the interesting property of leading to a prediction error which is linearly uncorrelated with the input \cite{bendat1993spectral,rocklin1985comparison}, therefore resulting in the best linear prediction of the output signal.

Once the transfer function $G_{IO}(f)$ is computed, its time-domain counterpart is recovered by means of an inverse Fourier transform:

\begin{equation}
\label{inversefourierfortheestimation}
g_{IO}(t)=\int_{-\infty}^{+\infty} G_{IO}(f) e^{-i 2 \pi f t} df.
\end{equation}

The function $g_{IO}(t)$ represents a linear convolution kernel, and therefore allows one to estimate the output via the convolution with the input signal

\begin{equation}
\label{convolutionforprediction}
O(t)=\int_{-\infty}^{\infty}g_{IO}(\tau)I(t-\tau)\rm d\tau.
\end{equation}

The dummy variable $\tau$ was introduced for the calculation of the convolution.

\subsection{Hillbert transform to evaluate causality}
\label{hilbert}

The convolution in equation \ref{convolutionforprediction} is taken between $\pm\infty$.
However, for a causal prediction, this operation should be limited from 0 to $+\infty$, so that only past information is used to predict the output signal.
The convolution operation is usually limited to the causal part of the convolution kernel, which permits an on-line prediction of the output.
This type of approach is necessary for active closed-loop applications, for example.
Nevertheless, the full convolution operation may also be used for data-reconstruction, a feature which was explored in \cite{sasaki2019transfer} for a turbulent boundary layer.

This implies that if $g_{IO}(t<0) \neq 0$, there is some reverse causality between input and output signals.
An evaluation of the convolution kernel in the time domain can therefore be used to determine the causality relation between input and output positions.
Such evaluations in the time domain could be cumbersome and, in some cases, qualitative, particularly for positions which are close to each other, where a feedback path usually occurs.
To avoid this problem, the evaluation of causality can be performed directly in the frequency domain by means of the Hilbert transform \citep{press2007numerical}.
This method has the advantage of testing several locations in a fast, computationally efficient manner.
The Hilbert transform is a linear operation that performs a convolution with the function $1/(\pi f)$,

\begin{equation}
\label{hilberneles}
H(F(f))=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{F(\Omega)}{f-\Omega}\mathrm{d}\Omega.
\end{equation}

Hilbert transforms are commonly used in signal processing for the design of causal filters, and the idea was adopted in flow control applications to determine suitable positions for the active control of shear layers \citep{sasaki2018_inpress}.
A causal filter results when $g_{IO}(t \leq 0)=0$.
This occurs in convectively unstable flows, such as a Tollmien-Schlichting or Kelvin-Helmholtz dominated shear flows, where upstream measurements are taken to estimate downstream fluctuations. Computation of the Hilbert transform, in this case, is made directly in the frequency-domain by means of a convolution with the previously calculated empirical transfer function.

The approach consists in exploring the following property of the transform: if the imaginary part of the frequency response (interpreted here as a transfer function) of a linear system is equal to the Hilbert transform of the real part, the system represents a causal, linear filter.
Hence, knowledge of the real part is sufficient to completely specify the system, with the imaginary part adding redundant information \cite{bendat2011random}.

Therefore, a quantitative evaluation of the causality of the transfer function $G_{IO}(f)$ may be obtained by computing the correlation between its imaginary part and the Hilbert transform of its real part. This parameter is given by

\begin{equation}
\label{pparamter}
P=\frac{\int_{-\infty}^{\infty}H(\Re[G_{IO}(f)])\Im[G_{IO}(f)]df}{\sqrt{\int_{-\infty}^{\infty}H(\Re[G_{IO}(f)])^2df}\sqrt{\int_{-\infty}^{\infty}\Im[G_{IO}(f)]^2df}}.
\end{equation}

Values of $P$ close to unity indicate a causal behaviour between input and output positions, with a reverse causality occurring as $P$ approaches zero.
This parameter will be used in the following section to evaluate the causality between two different locations.

\subsection{Spectral proper orthogonal decomposition}

In the current work, spectral proper orthogonal decomposition (SPOD) is applied in the streamwise/wall-normal plane.
As outlined in section III, the low-frequency shock dynamics can be considered mostly two-dimensional.
Therefore, spanwise averaging of the fourteen streamwise/wall-normal planes available from the LES data set is considered to isolate the $\beta=0$ fluctuations.
The effect of the spanwise averaging on the signal prediction framework discussed above is reported in section V (figure \ref{examplecorrelationsfasdf}).

The objective is to extract spatially-coherent fluctuations for a given frequency, in particular the low-frequency shock oscillations and its relation to the recirculating flow region.
The SPOD is used here as an auxiliary tool to observe the causes of the shock unsteadiness and to visualize the propagation of the fluctuations for a targeted frequency within the SBLI.
SPOD was employed in a number of studies \citep{picard2000pressure,cavalieri2013wavepackets,semeraro2016stochastic,towne2018spectral,sasaki2020role}, however, to the authors' knowledge, its application to the problem of a shock wave boundary layer interaction has not been reported previously.

The spectral decomposition is applied to the spanwise-averaged streamwise/wall-normal $(u,v)$ velocity fluctuations such that they are optimal modes to represent the turbulent kinetic energy. The modes are defined from the solution of the following integral equation:

\begin{equation}
\label{spodequation1}
\int \mathbf{\Gamma}(\mathbf{x},\mathbf{x'},f)\mathbf{W}(\mathbf{x'})\psi_j(\mathbf{x'},f)\mathrm{d}\mathbf{x'}=\lambda \psi_i(\mathbf{x},f).
\end{equation}

Here, $\mathbf{x}=(x,y)$ and $\mathbf{W}(\mathbf{x'})$ is a weighting matrix which allows one to account more either at arbitrary flow regions or due to the use of non-cartesian grid; $\psi_i$ corresponds to an eigenfunction ($i$th SPOD mode) with corresponding eigenvalue $\lambda$ and $\mathbf{\Gamma}(\mathbf{x},\mathbf{x'},f)$ is the two-point cross spectral density, which is defined from the Fourier transform of the correlation tensor

\begin{equation}
\label{definingthecrosspowerspectraltensor}
\mathbf{\Gamma}(\mathbf{x},\mathbf{x'},f)=\int_{-\infty}^{\infty}\mathbf{C}(\mathbf{x},\mathbf{x'},\tau)e^{i 2\pi f \tau}\mathrm{d}\tau.
\end{equation}

\noindent The correlation tensor $\mathbf{C}$ is obtained by:

\begin{equation}
\label{correlationtensor}
\mathbf{C}(\mathbf{x},\mathbf{x'},\tau)=\langle \mathbf{q}(\mathbf{x},t)\mathbf{q}^*(\mathbf{x}^\prime,t+\tau) \rangle,
\end{equation}

\noindent with $\mathbf{q}=(u,v)$, the velocity components, $*$ representing the conjugate transpose, and $\rangle . \langle$ supplying an ensemble time average, representing the expected value of a given realization of the flow field.

The method of snapshots was used to compute the SPOD modes, where the Welch's method is applied to compute the cross-spectral density.
A detailed description of this procedure can be found in \citet{towne2018spectral} or in the appendix of \citet{sasaki2020role}.
The discretized form of equation \eqref{spodequation1} is replaced by an eigenvalue problem \citep{towne2018spectral}, for each discrete frequency $f_k$, which reads:

\begin{equation}
\label{spodequation2}
\mathbf{\Gamma}_{f_k}\psi_{f_k}=\psi_{f_k}\lambda_{f_k},
\end{equation}

\noindent where $\mathbf{\Gamma}_{f_k}$ refers to the cross-spectral-density matrix and it is calculated via the ensemble averaging

\begin{equation}
\label{spodequation3}
\mathbf{\Gamma}_{ij}=\langle \hat{\mathbf{q}}_i \hat{\mathbf{q}}_j^* \rangle.
\end{equation}

The subscripts $f_k$ are not shown for simplicity, and $\hat{\mathbf{q_i}}=[\mathbf{u}_i, \mathbf{v}_i]^T$, with the superscript $T$ representing a transpose, denotes the $i$-th realisation of the velocity field (i.e., the $i$-th short-time Fourier transform of the velocity field in the Welch's method).

The elements in equation \ref{spodequation3} were determined by means of the Welch's method, with a triangular window and 75\% overlap of the segments.
Each segment has 5000 points with a time discretization of $\Delta t=0.026$.
The total number of elements in the time signal considered is of 342900, which results in 97 blocks for averaging.
These parameters were observed to be adequate to resolve the flow structures in the particularly low frequencies of interest.

\section{Transfer functions in the interaction zone}
\label{predictionresultsforthedata123}

\subsection{Prediction performance}

We start by assessing the accuracy of the single-input single-output linear transfer functions in predicting the velocity/pressure fluctuations in the average shock position.
The fluctuations were low-pass filtered using a finite-impulse-response (FIR) filter with a cut-off frequency of $St=0.3$ to avoid the influence of higher-frequency flow structures.
The filter is applied prior to the calculation of the transfer functions.
A window size consisting of 2000 points, which corresponds to $\Delta t=60$ with an overlap of 75\% was used in the Welch's method for calculation of the power spectra.
Similar results were also obtained for window sizes consisting of 5000 points.
To avoid spurious values of the transfer function, frequencies where $S_{IO}(f)$ is below a certain threshold are set to zero.
This occurs for frequencies higher than the cut-off filter.

Three input/output pairs were computed: pressure/pressure, velocity/pressure and velocity/velocity.
Here, only the causal part of the convolution kernel was considered to compute the prediction, \textit{i.e.} $g_{IO}$ was forced to zero when $t<0$.
To quantitatively evaluate the prediction performance, the normalized correlation between the estimation $O_{est}(t)$ and the LES data $O_{LES}(t)$ was computed

\begin{equation}
\label{equationforcorrelation}
C=\frac{\int_{-\infty}^{\infty}O_{LES}(t)O_{est}(t)\mathrm{d}t}{\sqrt{\int_{-\infty}^{\infty}O_{LES}^2(t)\mathrm{d}t}\sqrt{\int_{-\infty}^{\infty}O_{est}^2(t)\mathrm{d}t}},
\end{equation}

\noindent where $O_{LES}$ and $O_{est}$ are the LES and the estimated data, respectively.
The normalized correlation in the previous equation varies between zero and one, where one indicates a perfect correlation between predicted and LES signals.

Two input positions were considered separately, upstream and downstream of the shock.
The correlation was computed as a function of the output (prediction) streamwise position.
Figure \ref{examplecorrelationsfasdf} shows the resulting correlation between prediction and LES data, the vertical line highlights the input location.
The prediction was performed with and without the spanwise averaging of the input/output quantities.

\begin{figure}
\begin{center}
\subfigure[Input upstream of the shock, pressure/pressure estimation]{\includegraphics[width=0.48\textwidth]{./Figs/Updated_Figures/Prediction_Stream_Input_Fixed_p_Upshock-eps-converted-to.pdf}}
\subfigure[Input downstream of the shock, pressure/pressure estimation]{\includegraphics[width=0.48\textwidth]{./Figs/Updated_Figures/Prediction_Stream_Input_Fixed_p-eps-converted-to.pdf}}\\
\subfigure[Input upstream of the shock, velocity/pressure estimation]{\includegraphics[width=0.48\textwidth]{./Figs/Updated_Figures/Prediction_Stream_Input_Fixed_up_Upshock-eps-converted-to.pdf}}
\subfigure[Input downstream of the shock, velocity/pressure estimation]{\includegraphics[width=0.48\textwidth]{./Figs/Updated_Figures/Prediction_Stream_Input_Fixed_up-eps-converted-to.pdf}}\\
\subfigure[Input upstream of the shock, velocity/velocity estimation]{\includegraphics[width=0.48\textwidth]{./Figs/Updated_Figures/Prediction_Stream_Input_Fixed_Upshock-eps-converted-to.pdf}}
\subfigure[Input downstream of the shock,velocity/velocity estimation]{\includegraphics[width=0.48\textwidth]{./Figs/Updated_Figures/Prediction_Stream_Input_Fixed-eps-converted-to.pdf}}\\
\end{center}
   \caption{Correlation between prediction and LES for pressure and streamwise velocity components. The vertical dashed line presents the position of the input. Solid black lines correspond to the prediction performed with the spanwise averaging and dash-dot blue lines without it.}
  \label{examplecorrelationsfasdf}
\end{figure}

As expected, the correlations decay as the distance between input and output positions increases; for the cases where the same variable is used as the measurement and prediction, the correlations depart from large values close to unity.
When the input is upstream of the shock in figures 4(a,c,e), the correlation decays monotonically up to the mean shock location $x=0$ and stabilizes further downstream.
On the other hand, in the case of an input downstream of the shock wave shown in figures 4(b,d,f), the correlation also decays as the output position is moved further away from the input.
However, close to the shock location, there is a swift increase in the correlation values.
It should also be observed that even though there is a large separation between input and output, the correlation between prediction and input reaches values above 0.6.

The effect of the spanwise averaging is remarkable, leading to a significant increase in the correlation of the predicted signals, in particular when the input position is taken downstream of the shock motion.
This is in accordance with the spectra of figure \ref{2DSpectra}, which indicates a predominantly two-dimensional nature of the streamwise velocity shock-related fluctuations.
Hereafter, only spanwise averaged data is considered.

Figure \ref{examplepredictionstimedomain} presents a sample of the predicted streamwise velocity and pressure fluctuations at $x=0$ with input at $x_{in}=1.5$.
It is noticeable that, even though the amplitude of the LES signal is underestimated, both the phase and the low-frequency content of the signal (characteristic of the shock motion) are well captured by this approach, in spite of the large streamwise separation between input and output. The corresponding shape of the transfer function is depicted in figure \ref{exampleoftfdata}.

\begin{figure}
\begin{center}
\subfigure[Input downstream of the shock, pressure/pressure estimation]{\includegraphics[width=0.75\textwidth]{./Figs/Updated_Figures/Prediction_Comparison_Shock_pp-eps-converted-to.pdf}}\\
\subfigure[Input downstream of the shock, velocity/pressure estimation]{\includegraphics[width=0.75\textwidth]{./Figs/Updated_Figures/Prediction_Comparison_Shock_p-eps-converted-to.pdf}}\\
\subfigure[Input downstream of the shock, velocity/velocity estimation]{\includegraphics[width=0.75\textwidth]{./Figs/Updated_Figures/Prediction_Comparison_Shock_u-eps-converted-to.pdf}}\\
\end{center}
   \caption{Comparison between the time signals from LES and linear prediction, for pressure (a) and streamwise velocity (b) fluctuations. The predicted position corresponds to $x=0$ (average position of the shock) and the input is located at $x=1.5$ downstream of the shock.}
  \label{examplepredictionstimedomain}
\end{figure}

\subsection{Wall-normal behaviour and causality of the transfer function predictions}

To further investigate the prediction performance and to track the causal areas for the shock fluctuations, the output position was fixed at $x=0$, near the wall.
Transfer functions were computed considering a grid of positions along the streamwise and wall-normal directions as input signals.
As indicated above, we consider the spanwise averaged velocity fields low-pass filtered with the cut-off frequency $St=0.3$ prior to the computation of the transfer functions.
The same window size and overlap of the previous analysis were used for calculating the power spectra.
The pressure fluctuation was considered as the estimated variable, with the streamwise velocity fluctuations as the input.
Finally, the full convolution kernel including its non-causal counterpart was used here to evaluate signal reconstruction.

Figure \ref{pparameterandgroupvelocity}(a) presents the resulting correlations between the shock prediction and LES data for the different inputs considered.
High correlations appear in the near-shock region.
This is expected due to the coherent shock motion; such high correlations are due to an input related to shock oscillations predicting the fluctuations at the shock foot.
An interesting behavior occurs downstream of the shock ($x>1.0$), where correlations of the order of 0.8 start to occur, indicating that these can be used to accurately estimate the shock motion.
In the region $0.2<x<1.0$, within the recirculation bubble, the correlations drop to low values, a feature which can be assessed by the shape of the SPOD modes presented in following section.

In order to determine the causal effects of the shock, figures \ref{pparameterandgroupvelocity}(b,c) present two indices, the $P-$parameter, related to the causality between input and output positions, and the group velocity.
The latter can be estimated from the time delay observed in the empirically calculated transfer function and the streamwise separation between input and output positions.
As an example, consider the time-domain transfer function computed using the streamwise velocity fluctuation at the positions of $x_{in}=1.5$ and $x_{out}=0$, shown in figure \ref{exampleoftfdata}.
The streamwise separation is of -1.5, with a time delay for the peak of the transfer function corresponding to approximately $\Delta t=2$, which results in an estimated group velocity of $v_g=-0.75$, normalized with respect to the free-stream velocity $U_\infty$.
As the bulk of the transfer function is for positive time delays, the causality parameter $P$ is calculated as 0.9, indicating a nearly causal behavior.
Regions where the the $P$-parameter was between 0.3 and 0.7, or where the correlation $C$ between estimation and prediction was below 0.3, were not considered in the computation of the group velocity and the value was set to $v_g=0$.
This avoids a non-physical interpretation of this parameter near the $v_g=0$ line or where the transfer function estimation is inaccurate.

\begin{figure}
\centering
\includegraphics[width=0.75\textwidth]{./Figs/Updated_Figures/Example_TF-eps-converted-to.pdf}
\caption{Transfer function computed streamwise velocity fluctuation at the positions of $x_{in}=1.5$ and $x_{out}=0$.}
\label{exampleoftfdata}
\end{figure}

When both the correlation $C$ and $P-$parameter in figures \ref{pparameterandgroupvelocity}(a,b) are analyzed together, we note that the only causal input locations (large $P$) that lead to high correlations with shock motion (large $C$) surrounds the recirculation bubble and extends downstream of it.
In particular, the parameter $P$ and the group velocity indicate causality between shock oscillations at $x=0$ from fluctuations near the reattachment point and downstream of the recirculation bubble.
The computed negative group velocity points to fluctuations that are propagating upstream from the input location.
In addition, the high correlation values indicate that fluctuations at downstream locations precede shock oscillations, and enable an accurate prediction of its motion.

Fluctuations occurring upstream of the shock present a positive group velocity and are not strongly correlated with its motion.
These results are in agreement with previous observations from \citet{Touber2009}, who observed a negative phase speed in a similar configuration, and related it to the existence of a global mode.
The observations from \citet{Dupont2006} and \citet{Wu2008} also point to a strong correlation between low-frequency fluctuations and the reattachment region.
It is of interest to remark the existence of a region upstream of the shock with positive group velocity and large values of $P$.
This region contributes to causing the shock motion, which is in accordance with the experiments from \citet{Ganapathisubramani2007}. However, due to the relatively low correlation levels presented in figure \ref{examplecorrelationsfasdf}(a), this region has a relatively smaller contribution to the shock unsteadiness than the downstream zone in the present flow simulation.

\begin{figure}
\begin{center}
\subfigure[Correlations between estimated and LES data]{\includegraphics[width=0.80\textwidth]{./Figs/Updated_Figures/Correlations_OutShock_rev1-eps-converted-to.pdf}} \\
\subfigure[$P-$paramter]{\includegraphics[width=0.80\textwidth]{./Figs/Updated_Figures/P_parameter_OutShock_rev1-eps-converted-to.pdf}}\\
\subfigure[Estimated group velocity]{\includegraphics[width=0.80\textwidth]{./Figs/Updated_Figures/Group_Velocity_OutShocl_rev1-eps-converted-to.pdf}}
\end{center}
   \caption{Effect of the input position for prediction. Correlation values between predicted signal and LES data for the pressure fluctuations (a), $P-$parameter for computing the causal relationship (b) and streamwise velocity and estimated group velocity (c). The solid line indicates the recirculation bubble and the dashed line the boundary layer thickness. The cross indicates the considered output position}
  \label{pparameterandgroupvelocity}
\end{figure}

\section{Upstream travelling waves detected via SPOD}

\subsection{Shape and velocity of the spectral modes}

In order to better understand the spatial structure of the low-frequency oscillations, SPOD modes of the spanwise averaged $(u,v)$ velocity components were computed at $St=0.03$ and $St=0.3$.
The first Strouhal number corresponds to the low-frequency motion, and is representative of shock oscillations, whereas the second one is related to high-frequency boundary-layer and shear-layer fluctuations mostly uncorrelated to the shock.
Figure \ref{evaluationsofspodeigenvalues} depicts the resulting eigenvalues for these two selected frequencies.
In both cases, there is a clear dominance of the first mode with respect to the other spectral modes.
It is also noticeable that the low-frequency eigenvalues are higher over all the spectra.
Hence, we consider only the first mode in the following analysis.

\begin{figure}
\centering
\includegraphics[scale=0.65]{./Figs/Updated_Figures/Eigenvalues_SPOD-eps-converted-to.pdf}
\caption{Resulting SPOD eigenvalues computed for $St=0.03$ (blue circles) and $St=0.3$ (black squares).}
\label{evaluationsofspodeigenvalues}
\end{figure}

Figure \ref{SPOD_modes} presents the structure of the streamwise velocity fluctuations of the first SPOD mode for the frequencies $St=0.03$ and $St=0.3$.
The time evolution of the corresponding modes were obtained by multiplying the spectral mode by $e^{i 2 \pi f t}$.
The resulting snapshots are reported in the the supplemental movies.
Interestingly, a wave propagating with negative phase velocity (from right to left) is observed for $St=0.03$, on a position near the wall. Such motion surrounds the recirculation bubble toward the shock foot.
At a higher wall-normal position ($y \approx 0.15$) there is a wave generated on the shock, propagating with positive group velocity (from left to right).
The spatial footprint of this upstream travelling motion corresponds well to the causal locations presented in figure\ref{pparameterandgroupvelocity} (b).
This motion is schematized with the arrows in figure \ref{SPOD_modes}(a), and occurs for low-frequency modes in the vicinity of $St=0.03$.
For $St=0.3$, only a wave travelling with positive phase velocity appears.

\begin{figure}
\begin{center}
\includegraphics[scale=0.6]{./Figs/SPOD_modes.pdf}
\end{center}
   \caption{Snapshot of the first SPOD mode of the streamwise velocity component at $St=0.03$ (a) and $St=0.3$ (b). The arrows depict schematically the motion of the streamwise velocity fluctuations. The contours correspond to the real part of the first SPOD mode. The time evolution of both modes is reported in the supplemental material. }
  \label{SPOD_modes}
\end{figure}

To shed further light on the dynamics of the upstream and downstream travelling waves, an empirical dispersion relation is calculated.
This was obtained by computing the Fourier transform of the first SPOD mode ($u$) in the streamwise direction as a function of frequency $St$.
A similar approach was employed to study acoustic resonances in turbulent jets \citep{towne_cavalieri_jordan_colonius_schmidt_jaunet_bres_2017}.
The dispersion relation was computed for two wall-normal positions: $y=0.05$, close to the wall, but surrounding the recirculation bubble, and at $y=0.15$.
This procedure leads to the frequency-wavenumber power-spectral density diagram of figure \ref{empiricaldispersionrelationusingu}.
Finally, in order to better resolve the fluctuations occurring at different scales, the iso-contours are plotted in logarithmic scale.

\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.48\textwidth]{./Figs/Updated_Figures/SPOD_fftx_y005_rev1-eps-converted-to.pdf}}
\subfigure[]{\includegraphics[width=0.48\textwidth]{./Figs/Updated_Figures/SPOD_fftx_y015_rev1-eps-converted-to.pdf}} \\
\end{center}
   \caption{Empirical dispersion relations, computed from the first SPOD mode of the streamwise velocity at $y=0.05$ (a) and $y=0.15$ (b). The logarithmic of the resulting power spectral density was taken. The solid vertical lines indicate constant group velocities for the fluctuations.}
  \label{empiricaldispersionrelationusingu}
\end{figure}

The phase velocity of the fluctuations at a given $(St,k)$ pair is given by $v_p= 2\pi f/k$.
Therefore, waves with negative values of $k$ propagate upstream.
Since the group velocity reads $v_g=2\pi\partial f/\partial k$, we note that $v_g$ changes between the two considered wall-normal locations. For the low frequencies ($St<0.1$), there is a group of waves propagating upstream, with negative group velocity and this only occurs in the near-wall region.
On the other hand, at higher frequencies, the disturbances propagate with positive group velocity independently of the wall-normal position. The waves that travel with negative group velocity contribute to the causality relation between downstream areas and the shock motion, which was quantified by the correlations and the $P$-parameter obtained from the transfer functions.

\subsection{Reconstruction of the shock motion via SPOD}

In this section, we evaluate how much of the shock oscillations can be recovered from a limited set of data, specifically selected from regions upstream or downstream of the approximate shock position.
To quantify this analysis, we define the ratio

\begin{equation}
\label{equationforspodreconstruct}
\frac{|u_{SPOD_{cut}}|^2}{|u_{SPOD_{full}}|^2}=\frac{\int_{\Omega}\sum_{i=1}^{N_b}(\lambda_{cut_i}\psi_{cut_i})^2\mathrm{d}x}{\int_{\Omega}\sum_{i=1}^{N_b}(\lambda_{full_i}\psi_{full_i})^2\mathrm{d}x},
\end{equation}

\noindent where $\Omega$ is a horizontal line close to the wall encompassing the shock foot location, and the subscripts $full$ and $cut$ represent SPOD modes which were computed from the whole field and from a limited part of it using the weighting matrix $W$, respectively. Equation \ref{equationforspodreconstruct} represents the ratio between the energy of the streamwise velocity fluctuations from SPOD modes computed using partial knowledge of the flow field and all the data available.

When the whole field is used in the computation, the weighting matrix $W$ considers only the numerical quadrature corresponding to the use of a non-cartesian coordinate system.
To select specific flow regions, the numerical quadrature is multiplied by the matrix $W_1$, where $W_1=1$ for the selected region in the SPOD calculation, and $W_1=0$ in the region excluded from the calculation.
Therefore, the SPOD modes are computed considering only a portion of the field.
In spite of this fact, the whole field can be reconstructed; areas which were not accounted in $W_1$ but are observed in the SPOD modes can be related to fluctuations within the selected region where $W_1=1$.
This occurs due to the spatial coherence for the frequency where the SPOD mode is calculated.
We attempt to evaluate through SPOD in regions that exclude the shock itself but are coherent with the shock motion.
A similar approach was performed in \citet{sano2019trailing} to reconstruct the acoustic field scattered by an airfoil from knowledge of its near field fluctuations.

Figure \ref{exampleofareasinspodcalculation} provides distinct areas under evaluation.
The energy is integrated over $\Omega$, a straight line in the first grid position above the wall along $-0.05<x<0.05$.
This line is shown in the plot.
The spectral modes were calculated considering fluctuations within each of the rectangles limited by the dashed lines.
For each region, $W_1=1$ inside the rectangles, and zero outside them.
Two rectangles are highlighted, indicating regions I and II, a rectangle downstream (black) and upstream (blue) of the shock, respectively. In the following, these two regions are considered in detail.
The amount of energy that can be recovered is also shown in the figure \ref{exampleofareasinspodcalculation}.
It is clear that most of the energy of the shock fluctuations is related to downstream fluctuations.
When considering the upstream region only, less than 1\% of the energy is recovered.

\begin{figure}
\centering
\includegraphics[width=0.80\textwidth]{./Figs/Updated_Figures/SPOD_u_003_2_rev1_areas-eps-converted-to.pdf}
\caption{First SPOD mode for the streamwise velocity at $St=0.03$ and selected flow regions for the specific SPOD calculations. The portion of the energy at $\Omega$ which can be recovered in each case, in comparison to that of the full field, is highlighted in the corresponding rectangles.}
\label{exampleofareasinspodcalculation}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.7]{./Figs/Updated_Figures/Eigenvalues_Reconstruction_rev1-eps-converted-to.pdf}
\end{center}
   \caption{Eigenvalues of the SPOD calculation obtained from the full field and regions I/II at $St=0.03$.}
  \label{eigenvaluesoffre03and003reconstructed}
\end{figure}

The resulting eigenvalues of the SPOD modes at $St=0.03$ are shown in figure \ref{eigenvaluesoffre03and003reconstructed} for regions I and II.
The eigenvalues obtained from region I are close to the full data spectra, and approximately three orders of magnitude higher than those computed from region II, upstream of the shock.
The corresponding first SPOD modes (not shown here for brevity) show that shock wave fluctuations can be recovered from both upstream and downstream regions.
However, the amount of energy which is present in the shock is much lower when the upstream region is selected.
In addition, when the upstream area (region II) is considered, the spatial footprint of the upstream travelling wave becomes much thinner when compared to the downstream travelling wave.
These results provide further evidence of the relationship between the low-frequency motion, characteristic of the shock motion, and the flow regions downstream of it.

\section{Conclusions}

This study considered the computation of linear, single-input single-output transfer functions to detect flow regions responsible for the low-frequency shock oscillations in a SBLI interaction.
The transfer functions were used to predict unsteady shock fluctuations in a range of frequencies determined via low-pass data filtering ($St<0.3$).
It was observed that locations downstream of the shock could be used to accurately predict its fluctuations, where correlations of up to 0.6 were obtained between predicted and LES signals, even for the largest streamwise separation between input and output signals.
The correlation magnitude strongly depends on data spanwise averaging, pointing to a predominantly 2D mechanism governing the shock unsteadiness.

Evaluation of the transfer functions in the frequency domain by means of a Hilbert transform suggests a causal relationship between flow regions downstream of the shock and its motion.
When combined to the observed correlations, we conclude that the only causal inputs leading to high correlations with the shock are the locations surrounding and downstream of the recirculation bubble.
In these regions, the fluctuations exhibit a negative phase velocity.
These observations shed further light on the mechanisms suggested in previous investigations \cite{Touber2009,Piponniau2009}.
Causality was also detected between the upstream boundary layer and the shock motion, as reported in \citet{Ganapathisubramani2007}, but with a much lower correlation when compared to downstream locations.

The computation of the spectral modes from the SPOD analysis pointed to the existence of an upstream traveling wave in the leading mode at $St=0.03$.
The Fourier transform of the leading mode confirmed the negative group velocity for frequencies up to $St\approx 0.1$ in the near-wall region.
On the other hand, the same mode presented a positive group velocity at higher wall-normal positions.
In addition, we demonstrated that a significant portion of the energy of the shock motion can be recovered from a limited set of measurements taken downstream of the shock.
Finally, comparison of the corresponding eigenvalues indicated that the downstream region contributes more significantly to the shock unsteadiness than the upstream region.

\section{Acknowledgments}
This work was supported by the High Performance Computing resources of the Institut du Développement et des Ressources en Informatique Scientifique under the allocation 2014-2a1877 and 2015-2a1877 awarded by Grand Équipement National de Calcul Intensif for the LES computations. A.V.G.C. aknowledges funding from CNPq (grant 310523/2017-6).

\bibliographystyle{plainnat}