text:
\section{Introduction}
\label{intro}
Additive manufacturing (AM) is becoming increasingly popular in industries. Metal AM processes such as selective laser melting (SLM) \cite{YanSmith2016}, selective electron beam melting (SEBM) \cite{schwerdtfeger2014selective} and selective laser sintering (SLS) \cite{kruth2005binding} are layered approaches fabricating the products with complex shapes \cite{paul2014effect}. One of the bottlenecks for metal AM is the thermal stress caused by the high temperature gradients during the high-frequency heating/cooling cycles. High residual stress may result in defects \cite{li2018residual,taheri2017powder} such as the part distortion \cite{XIE2020102723}, delamination \cite{MUKHERJEE2017360} and poor fracture resistance \cite{MORI20066737,GrilliFracturePF2021,GRILLI2021100651}. The investigation of temperature history and stress evolution is therefore necessary for optimising the manufacturing parameters during the AM process. \\
Due to the large temperature gradients and high-frequency heating/cooling cycles during the manufacturing process, it is difficult to monitor the temperature and thermal stress in real time during experiments. For the residual stress measurement, the experimental approaches, such as the neutron diffraction \cite{grilli2020crystal,GrilliCOMPLAS2019}, X-ray diffraction \cite{HOCINE202030}, contour method \cite{vrancken2014residual} and semi-destructive hole drilling method \cite{mathar1934determination,martinez2019non}, are costly and sometimes limited to surface measurements. Numerical simulations are efficient to study the mechanical behaviour during the manufacturing process. The finite element method (FEM) has been widely used in simulating thermal stresses in AM \cite{ganeriwala2019evaluation}. One type of such models is based on the inherent strain method, where the deposited regions are activated with predefined strains to calculate the final deformation and residual stress of the part \cite{liang2018modified}. These models are computationally efficient, but the accuracy may be compromised due to the uncertainties of the inherent strains. Another type is thermo-mechanical coupling, where the temperature profile is calculated by solving the governing equations of heat transfer and then serves as the load for thermal stress calculation \cite{schoinochoritis2017simulation,QIAN2014312}.
Thermo-mechanical models are mostly constructed based on many geometrical simplifications and physical assumptions as the heat transfer part does not incorporate the thermal-fluid flow effects \cite{yan2020data}. Most of the FEM simulations failed to resolve the track morphology and molten pool evolution as all the tracks are assumed to be uniform. Therefore, those simulations are only able to predict the overall part distortion at the macroscopic scales, but cannot explain the thermal stress at the microscale. \\
An emerging type of the thermal stress models is the coupling between the computational fluid dynamics (CFD) and FEM models \cite{bailey2017laser,cheon2016thermal,chen2020high}. These coupled CFD-FEM models analyse the thermal stress and deformation using the high-accuracy temperature profiles output from the CFD simulations, particularly the powder-scale thermal-fluid flow simulations \cite{YAN2018210}, instead of those from simplified analytical models \cite{WEI2021100703}. Cheon et al. \cite{cheon2016thermal} proposed a CFD-FEM framework for welding to analyse the thermal stress along a single welding line, but the model did not consider the residual stress reset due to melting. Bailey et al. \cite{bailey2017laser} implemented the temperature and surface profiles from thermal-fluid flow simulations into an FEM residual stress model to simulate the multi-track manufacturing cases of the laser directed energy deposition (DED) process. The element activation and deactivation method in ABAQUS was used to simulate the melting and solidification, which is troublesome to input command lines to activate elements of the solidified materials in the input file and frequently encounters problems on numerical convergence and numerically-induced deformation of the newly-activated elements. In our previous work \cite{chen2020high}, we simulated the multi-layer multi-track SLM process by mapping the temperature profiles from the high-fidelity thermal-fluid flow simulation into the FEM model, enabling the high resolution of the molten pool evolution and the morphology (rough surfaces and inner voids) of the FEM model. The quiet element method was used to simulate the melting and solidification, where the melted materials are assigned with null material properties. \\
There are two critical but common problems in the aforementioned thermo-mechanical models and CFD-FEM models:
\begin{enumerate}
    \item How to represent the melting and solidification while avoiding numerically-induced deformation and numerical divergence. In the previous models, two methods are commonly used: first, the element activation method \cite{stender2018thermal,montevecchi2016finite,yang2016finite,lindgren1999simulation,lindgren2001modelling} can activate elements for solidified materials but cannot deactivate elements again due to remelting, while the surface nodal temperature obtained by the interpolation between active and inactive elements \cite{michaleris2014modeling} can cause significant inaccuracy and divergence problems. Second, the quiet element method \cite{lindgren2001modelling,michaleris2014modeling,ales2018integrated}, which is essentially the same as the stiffness degradation method commonly used in fracture mechanics \cite{Miehe2010,borden2018,GrilliDuarte2018,DuarteGrilli2018,GrilliKoslowski2019}, can remove the effect of the liquid phase on the solidified region, however upon solidification the liquid phase becomes solid with initial deformation that is incompatible with the deformation of the surrounding solid region, which induces un-physical deformation and may also cause numerical divergence.
    \item None of the previous models has considered the detailed grain structures and the resultant temperature-dependent anisotropic mechanical properties at the microscale. Using simplified homogenized constitutive laws, the previous models cannot resolve the micro-scale residual stress \cite{chen2019microscale}, which can result in unique mechanical properties of the AM parts significantly differing from those by the conventional manufacturing \cite{yan2020MRL}.
\end{enumerate}
In this paper, a crystal plasticity finite element model (CPFEM) \cite{GRILLI2015424,GrilliDIC2020} is developed, incorporating the temperature profile from the thermal-fluid flow simulation \cite{chen2020high} and grain structures from the grain growth simulation \cite{YangMin2021,yan2018integrated} and experimental data \cite{chen2019microscale}, as to be described in Section \ref{sec:cpmodel}. Particularly, an element elimination and reactivation method is proposed to model the melting and solidification. More importantly, the algorithm can reinitialise the state variables, as tested in a simple case in Section \ref{sec:testcase}, such as the plastic deformation of the solidifying elements, which is critical to avoid the pre-existing plastic deformation just after solidification. Simulation cases of SLM of 316 stainless steel are conducted as described in Section \ref{sec:simulationslm}. The proposed element elimination and reactivation method and the commonly used stiffness degradation method are compared in terms of prediction accuracy, convergence and computational efficiency. Based on the simulation results, the formation of residual stresses in different grains and the correlation between plastic deformation, depth and residual stress are discussed in Section \ref{sec:discussion}.

\section{Material model}
\label{sec:cpmodel}

\subsection{Crystal plasticity framework}
\label{sec:crystalplastframework}

A crystal plasticity material model including thermal response is used for this work.
The deformation gradient $\boldsymbol{F}$ can be decomposed into elastic, thermal and plastic parts \cite{asaro1977strain,kalidindi1998incorporation,vujovsevic2002finite,GRILLI2020104061}:
\begin{equation}
\boldsymbol{F}=\boldsymbol{F}_e \boldsymbol{F}_{th} \boldsymbol{F}_p \ ,
\label{eqn:elastoplasticdecomposition}
\end{equation}
where $\boldsymbol{F}_e$ is the elastic part including the stretching and rigid body rotation, $\boldsymbol{F}_{th}$ is thermal part and $\boldsymbol{F}_p$ is the irreversible plastic deformation, which evolves according to \cite{roters2010overview,ROTERS2019420}:
\begin{equation}
    \boldsymbol{L}_p = \dot{\boldsymbol{F}}_p \boldsymbol{F}_p^{-1} = \sum_{\alpha=1}^n \dot{\gamma}^{\alpha}\boldsymbol{m}^{\alpha} \otimes\boldsymbol{n}^{\alpha} \ ,
\end{equation}
where $\dot{\gamma}^{\alpha} $ is the shear rate of slip system $\alpha$. The vectors $\boldsymbol{m}^{\alpha}$ and $\boldsymbol{n}^{\alpha}$ are unit vectors which describe the slip direction and normal of slip plane $\alpha$, respectively. $n$ is the number of slip systems in active state. \\
The plastic strain rate on slip system $\alpha$ of face-centered cubic metals can be calculated by a power-law equation \cite{rice1971inelastic,peirce1983material}:
\begin{equation}
\label{eqn:gammadot}
\dot{\gamma}^{\alpha}=\dot{\gamma}_0 \left | \frac{\tau^{\alpha}}{\tau_c^{\alpha}(T)} \right |^{\frac{1}{m}} \rm{sgn}(\tau^{\alpha}) \ ,
\end{equation}
where $\tau^{\alpha}$ is the resolved shear stress, $\tau_c^{\alpha}$ is the temperature dependent critical resolved shear stress (CRSS), $\dot{\gamma}_0$ and $m$ are constants, which indicate the reference strain rate and rate sensitivity of slip, respectively. \\
An exponential model calibrated by experiments can be used to calculate the CRSS change with temperature decrease of each slip system \cite{daymond2006elastoplastic,ORNLpaper2020}:
\begin{equation}
\tau_c^{\alpha}(T)= \left ( k_A+k_{B}\cdot{\rm exp} [-k_{C}(T-T_0)] \right ) \cdot\tau_c^{\alpha}(T_0) \ ,
\end{equation}
where $\tau_c^{\alpha}(T)$ is the CRSS of slip system $\alpha$ at temperature $T$, $T_0$ is the reference temperature, $k_A$, $k_B$ and $k_C$ are constants. For slip system $\alpha$, its hardening behaviour is influenced by other slip systems $\beta$ as \cite{kalidindi1992approximate}:

\begin{equation}
\dot{\tau}_c^{\alpha} = \sum_{\alpha=1}^{n}h_{\alpha\beta} \left | \dot{\gamma}^{\beta} \right | \ ,
\end{equation}
where $h_{\alpha\beta}$ is the hardening matrix given by \cite{Diehl2017,IRASTORZALANDA2016184}:
\begin{equation}
h_{\alpha\beta} = q_{\alpha\beta} \left [ h_0 \left ( 1-\frac{{\tau_c^{\alpha}}}{\tau_\textrm{sat}} \right )^{a} \right ] \ ,
\end{equation}
where $q_{\alpha\beta}$ is a hardening coefficient matrix for latent hardening, $h_0$ is an initial hardening term, $\tau_{sat}$ is the saturation slip resistance and $a$ is a constant.
The parameters of the crystal plasticity model are calibrated using the experimental data in \cite{chen2019microscale}. \\
For the thermo-mechanical response during the laser scan process, the Green-Lagrange strain tensor can be expressed by the equation composed of thermal and elastic deformation gradients as:
\begin{equation}
\boldsymbol{E}_e = \frac{1}{2} \left ( \boldsymbol{F}_\textrm{th}^T\boldsymbol{F}_{e}^T\boldsymbol{F}_{e}\boldsymbol{F}_\textrm{th}-\boldsymbol{I} \right ) \ ,
\end{equation}
where $\boldsymbol{I}$ is the identity matrix. To describe how the size of the substance changes with the temperature, the linear thermal expansion coefficient $\alpha_l$ of 316L SS is introduced \cite{callen1998thermodynamics}. The volumetric thermal expansion coefficient can be expressed by $\alpha_v=3\alpha_l$ \cite{Dubrovinsky2002}. Since the thermal expansion coefficient of 316L stainless steel varies slightly with temperature  \cite{yadroitsev2015evaluation,grimvall1999thermophysical}, and the temperature dependence is given by fitting the experimental data into a linear equation \cite{jiang2012using}:
\begin{equation}
\alpha_v=\alpha_0+\alpha_1(T-T_0) \ ,
\label{eqn:linearTexpansion}
\end{equation}
where $\alpha_0$ denotes the thermal expansion coefficient at room temperature and $\alpha_1$ is the linear coefficient of the temperature dependence of $\alpha_v$. The infinitesimal volumetric change $dV$ from ${T_0}$ to ${T}$ can be described by \cite{vujovsevic2002finite}:

\begin{equation}
\alpha_{v}\mathrm{d}T= \frac{\mathrm{d}V}{V} \ ,
\label{eqn:infinitesimalTexpansion}
\end{equation}
where $V$ is the volume of the region and $T_0$ is the reference temperature (room temperature). Substituting equation (\ref{eqn:linearTexpansion}) into equation (\ref{eqn:infinitesimalTexpansion}), the thermal eigenstrain tensor $\boldsymbol{\alpha}$ can be expressed by \cite{grilli2018effect}:
\begin{equation}
\boldsymbol{\alpha}=\frac{1}{2}\left[{\rm exp}\left(\frac{1}{3}\alpha_{1}(T-T_{0})^2+\frac{2}{3}\alpha_{0}(T-T_0) \right) - 1 \right] \boldsymbol{I} \ .
\end{equation}
The 2nd Piola-Kirchhoff stress can be calculated by:
\begin{equation}
\boldsymbol{S}=\mathbb{C}(\boldsymbol{E_e}-\boldsymbol{\alpha}) \ ,
\label{eqn:pkstress}
\end{equation}
where $\mathbb{C}$ is the rank four elasticity tensor of 316L stainless steel \cite{Horiuci316SS,GRILLI2018104}. The elasticity tensor of the solid is temperature dependent \cite{daymond2006elastoplastic}, which can be described by a linear function for the components:
\begin{equation}
\mathbb{C}_{ij}^{\rm{solid}}(T)=\mathbb{C}_{ij}^{\rm{solid}}(T_0) + \frac{\mathrm{d}\mathbb{C}_{ij}^{\rm{solid}}}{\mathrm{d}T}(T-T_0) \ ,
\label{eqn:cijtempdep}
\end{equation}
where $\mathbb{C}_{ij}^{\rm{solid}}(T_0)$ is the elasticity tensor component at room temperature, and the derivative term indicates the change of the components caused by temperature change. Phase change can be modelled by changing the stiffness tensor $\mathbb{C}$, as to be explained in section \ref{sec:residualstiffmethod}. \\
Since positive strain represents tension and negative strain represents compression, the same applies to the stress tensor components. In the following, compression will appear as a negative stress component while tension will appear as a positive stress component. \\
The temperature $T$ is obtained from thermal-fluid flow simulations \cite{chen2020high}. It should be mentioned that the temperature mapping from thermal-fluid to mechanical simulation is one-way, and the heat transfer analysis is merely conducted in the high-fidelity thermal-fluid flow simulation, from which the thermal information are extracted for the crystal plasticity calculation. Thus, deactivation of the melted parts will not affect the temperature distribution in the representative volume.

\subsection{Residual stiffness method}
\label{sec:residualstiffmethod}

During the additive manufacturing process, the phase transitions cause significant change in the mechanical properties, which should be taken into consideration.
Since in our simulation the temperature profiles are obtained by mapping and interpolating the thermal-fluid flow simulation results to the finite element model in time and space \cite{YanLin2018,chen2020high}, there are liquid phase and gas phase existing in the model, which represent the melted material and gas above molten pool, respectively. Therefore, a residual stiffness method based on the temperature field is applied here to prevent sudden changes in mechanical properties caused by the temperature dependence of stiffness and avoid the convergence problem of simulation.
The specific method is to set a temperature range between the melting point $T_m$ of 316L stainless steel and the gas temperature $T_g$ set in the thermal-fluid flow simulation. In this temperature range, the stiffness tensor of the solid material is $\mathbb{C}_{ij}^{\rm{solid}}(T)$, which conforms to equation (\ref{eqn:cijtempdep}). \\
The temperature from the thermal-fluid flow simulation is recorded every $45$ $\mu$s and these time points are called the ``CFD output'' or ``temperature output''. The CPFEM simulations have smaller time steps, therefore an interpolation is necessary. For instance, in the time interval $\left [ t_1, t_2 \right ]$ the temperature is interpolated linearly from $T_1$ to $T_2$, as shown in Figure \ref{fig:cfdsteps}. \\
A phase identification algorithm has been implemented: if the temperature is higher than the melting temperature $T_m$ or lower than the gas temperature $T_g$, it means that the material is in liquid or gas phase. This algorithm transforms the stiffness matrix $\mathbb{C}_{ij}$ linearly with time between those of phase $1$ at $t_1$ and phase $2$ at $t_2$ in Figure \ref{fig:cfdsteps}. The formulas to find the stiffness tensor are listed in Table.\ref{tab:cfdsteps}. A residual stiffness coefficient $q_r$ is introduced to describe the residual stiffness tensor of the material:
\begin{equation}
\mathbb{C}_{ij}^{\rm{residual}}=q_{r}\cdot\mathbb{C}_{ij}^{\rm{solid}}(T_0) \ .
\label{eqn:residualstiff}
\end{equation}
The stiffness tensor $\mathbb{C}_{ij}$ in Table \ref{tab:cfdsteps} is used in equation (\ref{eqn:pkstress}) and the interpolation avoids discontinuities in  the stiffness tensor that would create convergence problems. \\
\begin{figure}[!htb]
\centering
\includegraphics[width=0.5\textwidth]{TimeInterpolationPhase}
\caption{\label{fig:cfdsteps} Schematic of the phase identification algorithm to implement the residual stiffness method for liquid and gas phases.}
\end{figure}
\begin{table}[!htb]
\normalsize
\caption{Time evolution of the stiffness tensor depending on the phase state at time $t_1$ and $t_2$.}
\label{tab:cfdsteps}
\centering
\def1.5{1.5}
\begin{tabular}{|c|c|c|}
\hline
Phase 1 & Phase 2 & Stiffness \\
\hline
Solid & Solid & $\mathbb{C}_{ij} = \mathbb{C}_{ij}^\mathrm{solid} \left ( T \right )$ \\
\hline
Solid & Gas/Liquid & $\mathbb{C}_{ij} = \frac{(t_2-t)}{(t_2-t_1)}\mathbb{C}_{ij}^\mathrm{solid} \left ( T_1 \right )
+ \frac{(t-t_1)}{(t_2-t_1)}\mathbb{C}_{ij}^\mathrm{residual}$ \\
\hline
Gas/Liquid & Solid & $\mathbb{C}_{ij} = \frac{(t_2-t)}{(t_2-t_1)}\mathbb{C}_{ij}^\mathrm{residual}
+ \frac{(t-t_1)}{(t_2-t_1)}\mathbb{C}_{ij}^\mathrm{solid} \left ( T_2 \right )$ \\
\hline
Gas/Liquid & Gas/Liquid & $\mathbb{C}_{ij} = \mathbb{C}_{ij}^\mathrm{residual}$ \\
\hline
\end{tabular}
\end{table}
The equations above are implemented into the MOOSE finite element framework \cite{permann2020moose}, where a Newton-Raphson approach to find the Cauchy stress tensor $\boldsymbol{\sigma}$ is applied. To obtain the 2nd Piola-Kirchhoff stress in equation (\ref{eqn:pkstress}), a function $\boldsymbol{\psi}=\boldsymbol{S}-\mathbb{C}(\boldsymbol{E_e}-\boldsymbol{\alpha})$ is defined, the absolute value of which must be minimised. An iterative return mapping algorithm, in which the stress $\boldsymbol{S}$ is the variable is implemented. The stress increment during the return mapping algorithm can be written as \cite{adhikary2016robust,chockalingam2013crystal}:
\begin{equation}
\boldsymbol{S}_{i+1}=\boldsymbol{S}_{i}-\left( \frac{\partial\boldsymbol{\psi}}{\partial\boldsymbol{S}} \right)^{-1}\cdot\boldsymbol{\psi} \ .
\end{equation}
where $i$ and $i+1$ indicate two subsequent iterations. The Cauchy stress $\boldsymbol{\sigma}$ is obtained by \cite{lee1969elastic}:
\begin{equation}
\boldsymbol{\sigma}=\frac{1}{\rm{det}(\boldsymbol{F_e})}\cdot\boldsymbol{F_e}\cdot\boldsymbol{S}\cdot\boldsymbol{F_{e}^T} \ .
\end{equation}
The Cauchy stress is used by the finite element solver to find the displacement field $\mathbf{u}$ that leads to stress equilibrium. \\
The parameters used in the simulations and the elastic constants of 316L stainless steel in Voigt notation are listed in Table \ref{tab:modelparameters}.

\begin{table}[!htb]
\normalsize
\caption{Material and model parameters used in the simulations \cite{daymond2006elastoplastic,chen2019microscale,jiang2012using,irastorza2017effect,Clausen1998}}
\label{tab:modelparameters}
\centering
\begin{tabular}{ll}
\hline\noalign{\smallskip}
Reference plastic strain rate ($\dot{\gamma}_0$) & $1 \times  10^{-7}$ s$^{-1}$  \\
Plastic strain exponent constant ($m$) & 0.1 \\
Hardening matrix ($h_{\alpha\beta}$) & 3839 MPa  \\
Hardening exponent ($a$) & 2.5  \\
Reference temperature ($T_0$) & 303 K \\
Melting temperature ($T_m$) & 1648.15 K  \\
Gas temperature ($T_g$) & 298.1 K \\
Critical resolved shear stress at $T_0$ ($\tau_c^{\alpha}$) & 213 MPa  \\
Saturation slip resistance ($\tau_{sat}$) & 302 MPa  \\
Volumetric thermal expansion coefficient at $T_0$ ($\alpha_0$) & 44.73$\times 10^{-6}$ K$^{-1}$ \\
Derivative of volumetric thermal expansion coefficient ($\alpha_1$) & 0.01011$\times 10^{-6}$ K$^{-2}$ \\
Temperature dependence of CRSS ($k_A$) & 0.53 \\
Temperature dependence of CRSS ($k_B$) & 0.47 \\
Constant in temperature dependence eq.(4) of CRSS ($k_C$) & 0.008 \\
Elastic constants at $T=T_0$ ($\mathbb{C}_{11}$) & 204.6 GPa \\
Elastic constants at $T=T_0$ ($\mathbb{C}_{12}$) & 137.7 GPa \\
Elastic constants at $T=T_0$ ($\mathbb{C}_{44}$) & 126.2 GPa \\
Derivative of elastic constant ($\mathrm{d}\mathbb{C}_{11}/{\mathrm{d}T}$) & -90.33 MPa K$^{-1}$ \\
Derivative of elastic constant ($\mathrm{d}\mathbb{C}_{12}/{\mathrm{d}T}$) & -45.10 MPa K$^{-1}$ \\
Derivative of elastic constant ($\mathrm{d}\mathbb{C}_{44}/{\mathrm{d}T}$) & -51.78 MPa K$^{-1}$ \\
Residual stiffness coefficient ($q_r$) & 0.01 and 0.001 \\
Laser scan speed & 700 mm/s \\
Laser power & 150 W \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table}

\subsection{Implementation of the element elimination and reactivation method}
\label{sec:elemelimreactimpl}

To model the melting and solidification of material in the SLM process, an element elimination and reactivation approach is employed. This approach is based on a moving subdomain paradigm that is originally developed \cite{Yushu2021ElemElim}, which is based on MOOSE's subdomain restricted system \cite{permann2020moose}. Specifically, the physical domain is divided into an active subdomain ($\Omega_a$) and an inactive subdomain ($\Omega_i$). The thermo-mechanical finite element computation is defined and carried out in $\Omega_a$, while $\Omega_i$ only carries the geometric information. In other words, the elements in the inactive domain, or inactive elements, do not contribute to the equation system that calculates stress equilibrium. \\
The active and inactive subdomains are determined only by the prescribed temperature field. At a typical FEM time step, any active element $\mathcal{T} \in \Omega_a $ with an average temperature that is above the melting temperature (i.e., $T_{\text{avg}} > T_m$) or below the gas temperature (i.e., $T_{\text{avg}} < T_g$) is moved from the active subdomain to the inactive one, i.e.:
\begin{equation}
\Omega_a \rightarrow \Omega_a \backslash \mathcal{T}, \quad \Omega_i \rightarrow \Omega_i \cup \mathcal{T}.
\end{equation}
 This process mimics the solid material being eliminated from the physical domain, which happens when the material is melted around the laser beam, and is therefore referred to as element elimination.
Similarly, when the average temperature of one inactive element $\mathcal{T} \in \Omega_i $ falls to the solid temperature range (i.e., $T_g \leq T_{\text{avg}} \leq T_m$), it is moved to the active subdomain, i.e.:
\begin{equation}
\Omega_a \rightarrow \Omega_a \cup \mathcal{T}, \quad \Omega_i \rightarrow \Omega_i \backslash \mathcal{T}.
\end{equation}
This mimics the process of material being added back to the solid physical domain, which happens when the material solidifies behind the laser beam,  thus is referred to as element reactivation.

After elements being transferred between the two subdomains, several tasks need to be accomplished before the subsequent computation. First, the boundary information is updated based on the coverage of the new subdomains. Here, the subdomain boundaries are updated for both $\Omega_a $ and $\Omega_i$. The boundary conditions are projected to the updated boundaries correspondingly.
Second, the initial conditions of the displacement variables and the material properties that define the crystal plasticity behaviour, such as the plastic deformation (see Section \ref{sec:crystalplastframework}), are projected to the quadrature points of the reactivated elements. Projection of the initial condition indicates that the displacement, strain, and stress calculations are reset to zero. This allows for the restart of the plastic deformation calculation during the simulation, which cannot be implemented easily in the stiffness degradation method \cite{Duarte2019APS}. Indeed, if the elements are not eliminated and reactivated, they can still carry a certain amount of load. In the stiffness degradation method, a sudden change of the plastic deformation due to melting and solidification causes convergence problems because an elastic deformation jump, and consequent stress discontinuity, occur. \\
Finally, the element elimination and reactivation method resolves the issue that is brought by large deformation from the reactivated elements if they are set to carry the displacement history upon reactivation.

\section{Test case: element elimination during tension}
\label{sec:testcase}

The element elimination and reactivation method is tested on a simple example geometry, as shown in Figure \ref{fig:methodtest}. The aim of this test is to demonstrate that inactive elements do not contribute to the stress equilibrium calculations, i.e. they do not exert stress on the active elements, and to prove that boundary conditions are implemented correctly on the continuously changing free surface. \\
A cubic representative volume is constituted of 10 $\times$ 10 $\times$ 10 elements and has a size of 1 $\mu$m $\times$ 1 $\mu$m $\times$ 1 $\mu$m. Boundary condition $\mathbf{u}_x = 0$ is applied on the surface $x = 0$, $\mathbf{u}_y = 0$ on $y=0$ and $y = 1$ $\mu$m, $\mathbf{u}_z = 0$ on $z=0$ and $z = 1$ $\mu$m. Displacement boundary condition is imposed on the surface $x = 1$ $\mu$m. $\mathbf{u}_x$ is increased up to $0.001$ $\mu$m on that surface before the element elimination algorithm starts. \\
The element elimination algorithm is based on a prescribed temperature field, which is chosen in such a way that elements from the surface $z = 1.0$ $\mu$m are progressively removed, down to the surface $z = 0.5$ $\mu$m. The displacement $\mathbf{u}_x$ on the surface $x = 1$ $\mu$m is kept constant during element elimination, as shown in Figure \ref{fig:methodtest}(a)-(c), therefore $\boldsymbol{\varepsilon}_{xx}$ remains constant at the value of $0.001$. \\
Plasticity and thermal eigenstrain are not activated in this simulation case, and the elastic constants are set as $\mathbb{C}_{11} = 1.5$ GPa and $\mathbb{C}_{12} = 0.75$ GPa. $\boldsymbol{\varepsilon}_{yy}$ and the off diagonal strain components are zero. When no element elimination has taken place, the boundary condition $\mathbf{u}_z = 0$ on the surface $z=1$ $\mu$m holds, therefore $\boldsymbol{\varepsilon}_{zz} = 0$ and $\mathbf{u}_z = 0$ everywhere in the representative volume, as shown in Figure \ref{fig:methodtest}(d). This lateral constraint is reflected by the large positive value of the stress component $\boldsymbol{\sigma}_{zz}$ in Figure \ref{fig:methodtest}(g). Therefore, a tensile stress is present on that surface. \\
\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[width=0.35\textwidth]{DispXTime5}}
\subfloat[(b)]{\includegraphics[width=0.35\textwidth]{DispXTime20}}
\subfloat[(c)]{\includegraphics[width=0.35\textwidth]{DispXTime90}}
\newline
\subfloat[(d)]{\includegraphics[width=0.38\textwidth]{DispZTime5}}
\subfloat[(e)]{\includegraphics[width=0.38\textwidth]{DispZTime20}}
\subfloat[(f)]{\includegraphics[width=0.38\textwidth]{DispZTime90}}
\newline
\subfloat[(g)]{\includegraphics[width=0.38\textwidth]{SigmaZZTime5}}
\subfloat[(h)]{\includegraphics[width=0.38\textwidth]{SigmaZZTime20}}
\subfloat[(i)]{\includegraphics[width=0.38\textwidth]{SigmaZZTime90}}
\caption{\label{fig:methodtest} Test case for the element elimination and reactivation method. (a),(b),(c) Firstly, a displacement $\mathbf{u}_x$ is applied and elements are eliminated. (d) The lateral displacement starts at zero, (e)-(f) then lateral contraction appears due to the newly created free surface perpendicular to the $z$ axis. For the same reason, (g) the lateral stress $\boldsymbol{\sigma}_{zz}$ starts at 750 MPa, (h)-(i) then $\boldsymbol{\sigma}_{zz}$ becomes zero.}
\end{figure}
After the element elimination algorithm starts, the top surface perpendicular to the $z$ axis becomes unconstrained. Therefore, lateral contraction takes place because of the constant displacement along the $x$ axis. The lateral contraction can be calculated as follows: $\boldsymbol{\varepsilon}_{zz} = - \left ( \mathbb{C}_{12} / \mathbb{C}_{11} \right ) \boldsymbol{\varepsilon}_{xx} = -0.0005$. For example, in Figure \ref{fig:methodtest}(e), 8 elements are still present along the $z$ axis, therefore the lateral displacement on the top surface is given by $\mathbf{u}_z = \boldsymbol{\varepsilon}_{zz} \cdot 0.8$ $\mu$m $= -0.0004$ $\mu$m. The same holds in Figure \ref{fig:methodtest}(f), in which the lateral displacement on the top surface is $\mathbf{u}_z = -0.00025$ $\mu$m because 5 elements are still present along the $z$ axis. \\
As shown in Figures \ref{fig:methodtest} (h)-(i), there is no lateral stress $\boldsymbol{\sigma}_{zz}$ after the element elimination algorithm starts. This proves that the method developed creates new free surfaces in the representative volume and the eliminated elements do not apply stress on the active elements.

\section{Simulation of selective laser melting}
\label{sec:simulationslm}

\subsection{Simulation setup}
\label{sec:simulationsetup}

Crystal plasticity simulations of selective laser melting are carried out using the representative volume in Figure \ref{fig:grains}. A single scan track is modelled. The dimensions along the X, Y and Z axes are $200$ $\mu$m, $160$ $\mu$m and $108$ $\mu$m respectively. The FEM model represents a small region inside a larger sample. However, it is representative of the texture because of the large number of grains included. \\
The temperature field is obtained from thermal-fluid flow simulations
\cite{yan2017multi,wang2020evaporation,chen2020high}. The governing equations are the conservation of mass, momentum and energy. Most of the physical factors are incorporated, including the ray-tracing heat source model for laser, thermal conduction, surface radiation and convection, latent heat, evaporation, recoil pressure, viscosity, buoyancy force, surface tension and Marangoni effect. The details of the model and material parameters are reported in our previous works \cite{yan2017multi,wang2020evaporation}.  \\
The grain structure is obtained from a phase field model for grain growth \cite{YangMin2021}, where heterogeneous nucleation \cite{simmons2000phase} and the initial grain structures in powder particles and substrate are incorporated, and the same temperature profile from the thermal-fluid flow simulation is used. The phenomena, including grain nucleation and growth, competitive growth, epitaxial growth from powder particles and substrate, and grain coarsening in heat affected zones (HAZs), are comprehensively considered and experimentally validated. The details of the model and parameters are reported in our previous work \cite{YangMin2021}. \\
A total of 4522 grains are present with different orientations. The voxels obtained from the phase field model for grain growth are mapped directly into the FEM mesh, therefore the shape of the grains is reproduced accurately. The Euler angles $\varphi_1$, $\Phi$ and $\varphi_2$ are assigned based on electron backscatter diffraction (EBSD) experiments provided by Dr. Yin Zhang and published elsewhere \cite{chen2019microscale} to reproduce the realistic texture of AM 316 stainless steel. This is important because the grain orientation determines both the macroscopic and microscopic mechanical properties of the sample. On the other hand, the 3D shape of the grains cannot be fully extracted from EBSD measurements, therefore simulation results from the grain growth model are adopted. Figure \ref{fig:grains} shows the value of $\varphi_1$. The average grain size is approximately $16$ $\mu$m. However, grains in the centre of the laser track tend to be elongated along the $z$ axis, as shown in Figure \ref{fig:grains}. \\
The temperature field obtained from thermal fluid-flow simulations is projected on the mesh of the FEM simulations. Melting in the crystal plasticity simulations is described either by the stiffness degradation model described in section \ref{sec:residualstiffmethod} or by the element elimination and reactivation method described in section \ref{sec:elemelimreactimpl}. The temperature fields at different time steps are shown in Figure \ref{fig:temperature} (a)-(c), in which the eliminated elements are hidden from view. A structured mesh made of first order hexahedral elements is used. Simulations are carried out with the element size being $4$ $\mu$m or $8$ $\mu$m, respectively. The numbers of elements are 54000 and 7000 respectively. \\
\begin{figure}[!htb]
\centering
\includegraphics[width=0.6\textwidth]{GrainsLength3.png}
\caption{\label{fig:grains} Representative volume for the selective laser melting simulations and grain structure implemented from the phase field simulation for the crystal plasticity model.}
\end{figure}
The boundary conditions used in the simulations represent the constraints that are applied by the substrate and by the surrounding material. The displacement components are zero on the surfaces $x = 0$, $x = 200$ $\mu$m, $y = 0$ and $y = 160$ $\mu$m. The vertical displacement $\mathbf{u}_z$ is zero on the bottom surface, $z = 0$, to avoid translation along the $z$ axis. The displacement is free on the upper surface ($z = 108$ $\mu$m), where the laser beam is present.

\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[height=0.3\textwidth]{Temperature3NoScale}}
\subfloat[(b)]{\includegraphics[height=0.3\textwidth]{Temperature8NoScale}}
\subfloat[(c)]{\includegraphics[height=0.3\textwidth]{Temperature21Scale}}
\newline
\subfloat[(d)]{\includegraphics[height=0.3\textwidth]{SigmaXXtime99}}
\subfloat[(e)]{\includegraphics[height=0.3\textwidth]{SigmaXXtime226}}
\subfloat[(f)]{\includegraphics[height=0.3\textwidth]{SigmaXXtime451}}
\caption{\label{fig:temperature} Snapshots of the temperature field using the element elimination and reactivation method at (a) $t = 99$ $\mu$s, (b) $t = 226$ $\mu$s, (c) $t = 451$ $\mu$s; (d)-(f) corresponding stress component $\boldsymbol{\sigma}_{xx}$ using the stiffness degradation method.}
\end{figure}

\subsection{Simulation convergence}
\label{sec:simulationconvergence}

The convergence speed of the simulations using the stiffness degradation method depends mainly on the residual stiffness parameter $q_r$. A value smaller than 0.001 significantly increases the simulation time, but there is very small difference in the plastic deformation between the simulations with $q_r=0.01$ and $q_r=0.001$, as will be shown in the following.
Simulations with both $q_r=0.01$ and $q_r=0.001$ are carried out. \\
If the element elimination method is used, the time step must be decreased immediately before and after the temperature steps at which element elimination takes place, as described in section \ref{sec:cpmodel}. In this way, good convergence can be obtained. Specifically, the time step is reduced to $0.1$ $\mu$s in an interval of $2$ $\mu$s around each temperature step, i.e. every $45$ $\mu$s. For instance, considering the first temperature step, the time step is reduced to $0.1$ $\mu$s at $t = 44$ $\mu$s and it is increased again to $1$ $\mu$s after $t = 46$ $\mu$s.
This time step reduction applies only during the laser scan, while time step is kept constant at $1$ $\mu$s during the subsequent cooling stage, which is the second part of the simulation, carried out after the laser beam has exited the geometry. By contrast, in the stiffness degradation method, the value of the time step is always $1$ $\mu$s.
Since the eliminated elements do not contribute to the residual or Jacobian calculation of the solver, the computation of each time step is slightly faster if the element elimination method is used. However, a reduction of the time step size before and after the element elimination steps is necessary in that case. Therefore, the simulations carried out with the element elimination method are not necessarily faster, unless the number of eliminated elements is a significant fraction of the total number of elements. This is shown by the computational costs reported in table \ref{tab:simulationtime}. \\
\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[height=0.38\textwidth]{UzTime771ElemDel}}
\subfloat[(b)]{\includegraphics[height=0.28\textwidth]{ElementMemory3D}}
\newline
\subfloat[(c)]{\includegraphics[height=0.38\textwidth]{UzTime765_2}}
\subfloat[(d)]{\includegraphics[height=0.38\textwidth]{UzTime771}}
\caption{\label{fig:convergence} Vertical displacement $\mathbf{u}_z$ at $t = 771$ $\mu$s: (a) just before element reactivation, (c)-(d) just after element reactivation in case the reinitialisation $\mathbf{u}_z = 0$ or $\mathbf{u}_z = 3$ $\mu$m is used respectively. (b) Element distortion due to the reinitialisation $\mathbf{u}_z = 0$.}
\label{fig:convergenceuz}
\end{figure}
\begin{table}[!htb]
\normalsize
\caption{Computational costs of the simulations.}
\label{tab:simulationtime}
\centering
\begin{tabular}{|c|c|c|c|}
 \hline
 Mesh size & CPUs & Element elimination & Stiffness degradation \\
 \hline
 4 $\mu$m & 16 & 99 hours & 79 hours \\
 \hline
 8 $\mu$m & 12 & 9.5 hours & 9.5 hours \\ \hline
\end{tabular}
\end{table}
As stated before, the element elimination steps require a reduction of the time step, however, they do not cause further convergence problems. The time steps at which element reactivation takes place leads to the following problem. Because of thermal expansion, there is a significant displacement along the $z$ axis on the top surface, as shown in Figure \ref{fig:convergenceuz}(a). Normally, reactivated elements have zero displacement because the displacement variable is reinitialised. In the MOOSE framework, the displacement in the reinitialised nodes is assigned in such a way as to get an average zero displacement at the eight integration points, as shown in Figure \ref{fig:convergenceuz}(b). Therefore, when an element is reactivated with zero displacement, the excessive distortion leads to divergence, as shown in Figure \ref{fig:convergenceuz}(c). \\
Several attempts have been carried out in which the initial conditions on the displacement after element reactivation are changed as a function of the coordinates. Surprisingly, a simple solution is to reinitialise $\mathbf{u}_z = 3$ $\mu$m when elements are reactivated, while $\mathbf{u}_x$ and $\mathbf{u}_y$ are reinitialised to zero everywhere. This displacement value is closer to the average displacement of the molten pool region on the top surface in the simulation using the residual stiffness method, while it also leads to good convergence. Other reinitialisation values for $\mathbf{u}_z$, such as 1 $\mu$m, 2 $\mu$m, 5 $\mu$m and $6$ $\mu$m have been used, but they have lead to excessive element distortion just after reinitialisation. Only values between 3 $\mu$m and 4 $\mu$m remove the element distortion problem, as shown in Figure \ref{fig:convergenceuz}(d). This implementation is made possible by the MOOSE framework because the initial conditions can be set as a function of time. At $t = 0$, the displacement is initialised to zero, while for $t > 0$, the components $\mathbf{u}_x$, $\mathbf{u}_y$ are reinitialised to zero and $\mathbf{u}_z$ is reinitialised to $3$ $\mu$m everywhere in the geometry. This strategy can be applied to any additive manufacturing simulation and can be extended by setting the initial conditions as a more complex function of space and time.

\subsection{Comparison between the computational methods}
\label{sec:simulationresults}

The two computational methods presented in sections \ref{sec:residualstiffmethod} and \ref{sec:elemelimreactimpl} are compared. The plastic deformation is a very important quantity because it determines the residual stress after the laser scan. The diagonal components of $\boldsymbol{F}_p-\boldsymbol{I}$ are averaged over the central part of the representative volume, constituted of the highlighted elements in Figure \ref{fig:plasticstraincentre}(d). This procedure is carried out for simulations using the element elimination method and the stiffness degradation method. \\
Figures \ref{fig:plasticstraincentre}(a)-(c) show the comparison of simulation results by the two methods. During laser scan, the hot region undergoes thermal expansion and applies compression on the neighbouring regions. Therefore, plastic compression along the $x$ and $y$ axes takes place in the central region. Since plastic deformation is isochoric, the elements expand along the $z$ axis, as shown in Figure \ref{fig:plasticstraincentre}(c). \\
The plastic deformation is always larger using the element elimination and reactivation method. If the stiffness degradation method is used, the residual plastic deformation component along the $y$ axis, after laser scan and cooling, is underestimated by about 20\%. The difference is particularly important at time $400$-$500$ $\mu$s, when the laser just passes the representative volume but the melting pool is still deep. In this time interval, the stiffness degradation model underestimates the plastic deformation by a factor 2. \\
\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[height=0.4\textwidth]{FpxxCentre}}
\subfloat[(b)]{\includegraphics[height=0.4\textwidth]{FpyyCentre}}
\newline
\subfloat[(c)]{\includegraphics[height=0.4\textwidth]{FpzzCentre}}
\subfloat[(d)]{\includegraphics[height=0.32\textwidth]{SelectionCentre}}
\caption{\label{fig:plasticstraincentre} (a)-(c) Plastic strain components in the centre of the representative volume, averaged over the highlighted elements in (d).}
\end{figure}
Overall, in the stiffness degradation method, the presence of low stiffness elements prevents the elements near the melting pool to accommodate part of the plastic deformation. Simulations with both $q_r=0.01$ and $q_r=0.001$ are carried out and the plastic deformation along the $y$ axis is shown in Figure \ref{fig:Fpqr} when the laser beam is approximately in the centre of the representative volume. The strong similarity between the plastic deformation in the two simulations confirms that the range of values chosen for the parameter $q_r$ does not affect the simulation results. The other components of the plastic deformation gradient are also very similar, therefore a value $q_r = 0.001$ is suitable for this kind of simulations.

\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[height=0.41\textwidth]{FpyyTime171usgg0p01.png}}
\subfloat[(b)]{\includegraphics[height=0.41\textwidth]{FpyyTime170usgd0p001.png}}
\caption{\label{fig:Fpqr} Component $\left ( \boldsymbol{F}_p \right )_{yy} - 1$ of the plastic deformation gradient at time $t = 179$ $\mu$s using the residual stiffness method with (a) $q_r=0.01$ and (b) $q_r=0.001$.}
\end{figure}

\subsection{Grain orientation and plastic deformation}
\label{sec:grainoripdef}

The element elimination and reactivation method will be used in the following simulation results to understand the correlation between plastic deformation, residual stress and grain orientation. Firstly, the correlation between residual plastic deformation (after the laser scan and cooling) and grain orientation is analysed. The scatter plots in Figure \ref{fig:schmidfactor}(a)-(c) are obtained as follows: the plastic deformation components are averaged over each element in the representative volume and one element is selected every 16 $\mu$m and represented as a point in the scatter plots. The distribution of these elements is uniform in the representative volume and this algorithm allows to select a representative set of elements that do not fill completely the scatter plot. Given the Euler angles in each selected element, the maximum Schmid factor among the slip systems is calculated for load along a particular direction ($x$, $y$ or $z$ axes respectively) \cite{Grilli2016thesis}. This is motivated by the fact that, despite of the complex load undergone by different regions, the stress component along one direction should in principle induce plastic deformation along that particular direction. The colour of the points in Figure \ref{fig:schmidfactor}(a)-(c) corresponds to the position of the element along the $z$ axis. \\
\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[height=0.5\textwidth]{FpxxSchmidFactorCoordsZLessP}}
\subfloat[(b)]{\includegraphics[height=0.5\textwidth]{FpyySchmidFactorCoordsZLessP}}
\newline
\subfloat[(c)]{\includegraphics[height=0.5\textwidth]{FpzzSchmidFactorCoordsZLessP}}
\caption{\label{fig:schmidfactor} Correlation between maximum Schmid factor among the slip systems and plastic deformation after laser scan.}
\end{figure}
\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[height=0.42\textwidth]{FpxxRangeSchmidFactor}}
\subfloat[(b)]{\includegraphics[height=0.42\textwidth]{FpyyRangeSchmidFactor}}
\newline
\subfloat[(c)]{\includegraphics[height=0.42\textwidth]{FpzzRangeSchmidFactor}}
\caption{\label{fig:rangeschmidfactor} Correlation between maximum Schmid factor among the slip systems and the range of plastic deformation after laser scan.}
\end{figure}
As shown in Figure \ref{fig:schmidfactor}, the maximum plastic deformations occur on some points that are at an intermediate depth (about 40 $\mu m$ to 80 $\mu m$). This is related to the solidification and thermal stress around the molten pool, as also revealed in Figure \ref{fig:FpMeltPool}. Moreover, it is clear that a higher number of grains with larger maximum Schmid factors along a specific direction ($x$, $y$ or $z$) tends to have larger plastic deformation along that direction. However, there are grains with large maximum Schmid factors that show little plastic deformation. Many of this kind of points are dark blue in Figure \ref{fig:schmidfactor}(a)-(c), which means they are near the substrate. Because the substrate is fixed during simulation, the plastic deformation near the substrate is limited. Nevertheless, there are still some points at intermediate depth that do not have much plastic deformation, this is quite surprising and reflects the complexity of the thermal load in different regions. This behaviour will be clarified in the following, when the plastic deformation during laser scan is shown near the melting pool. There are also some grains with larger Schmid factor that show little plastic deformation and are located near the top stress-free surface (dark red points), as shown in Figure \ref{fig:schmidfactor}(a)-(c). \\
To investigate the directionality of the plastic deformation, the range of the plastic deformation components is investigated, which shows a greater correlation with the maximum Schmid factor \cite{IRASTORZALANDA2017157}, as shown in Figure \ref{fig:rangeschmidfactor}(a)-(c). The range is calculated as the difference between maximum and minimum values of the plastic strain components for a given maximum Schmid factor.
This correlation does not depend strongly on the depth $z$.
The majority of the grains undergo plastic compression along the $x$ and $y$ axes, while plastic tension is observed along the $z$ axis. This is consistent with the average plastic strain in Figure \ref{fig:plasticstraincentre}.

\subsection{Plastic deformation and residual stress}
\label{sec:pdefresstress}

The residual stress, after laser scan and cooling, is due to the heterogeneity of the plastic strain in different grains. In order to maintain strain compatibility, the elastic deformation must compensate the difference in plastic deformation, therefore residual stresses are generated, as modelled by equations (\ref{eqn:elastoplasticdecomposition}) and (\ref{eqn:pkstress}).
The correlation between residual stress and plastic deformation
can be visualised using the scatter plots in Figure \ref{fig:resstress}.
Each point represents an element in the representative volume, selected every 16 $\mu$m,
with its corresponding residual stress and plastic strain, averaged over the integration points.
The colour of each point represents the $z$ coordinate of the corresponding element. \\
It is certainly true that there is a correlation between residual stress and plastic strain.
Compressive plastic strain along the $x$ and $y$ axes corresponds to tensile residual stress components
along those directions, while tensile plastic strain along the $z$ axis corresponds to
compressive residual stress along $z$. However, few points show the opposite behaviour: for instance, a tensile residual stress $\boldsymbol{\sigma}_{yy}$ is associated with a small tensile plastic strain
along the $y$ axis in some regions. Overall, the tensile residual stress components $\boldsymbol{\sigma}_{xx}$ and $\boldsymbol{\sigma}_{yy}$ are similar and large, therefore these pre-tensions can induce fracture along both the $x$ and $y$ directions. \\
\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[height=0.5\textwidth]{ResStressXXCoordsZLessP.png}}
\subfloat[(b)]{\includegraphics[height=0.5\textwidth]{ResStressYYCoordsZLessP.png}}
\newline
\subfloat[(c)]{\includegraphics[height=0.5\textwidth]{ResStressZZCoordsZLessP.png}}
\caption{\label{fig:resstress} Correlation between the components of the residual stress and components of the plastic deformation gradient.}
\end{figure}
The plastic deformation depends also on the depth. In Figure \ref{fig:resstress}, it is evident that
the regions at the top of the representative volume, represented by dark red points,
are not the ones accommodating the largest plastic strain. The points with the largest
plastic deformation are the ones with depth in the interval 30 $\mu$m $< z < 80$ $\mu$m.
These are the points that are closer to the bottom of the melting pool during laser scan.
However, the points at the top surface can have residual stresses that are similar
to points that undergo larger plastic deformations. This shows that the influence of the plastic deformation in the substrate is very important; it will be discussed in more details in section \ref{sec:discussion}. It is important to note that some elements show values of the residual stress components that are larger than the critical resolved shear stress to induce plastic deformation. This is due to large volumetric stress in those elements. Moreover, some dark red points in Figure \ref{fig:resstress}(c) have residual stress $\sigma_{zz}$ that is different from zero. These points are picked just below the surface, in the coordinate interval 90 $\mu$m $< z < 100$ $\mu$m. Indeed the stress component $\sigma_{zz}$ on the free top surface is very close to zero, as shown in Figure \ref{fig:SigmaFpMeltPool}(c). \\
\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[height=0.41\textwidth]{FpxxTime179.png}}
\subfloat[(b)]{\includegraphics[height=0.41\textwidth]{FpyyTime179.png}}
\newline
\subfloat[(c)]{\includegraphics[height=0.41\textwidth]{FpzzTime179.png}}
\caption{\label{fig:FpMeltPool} Components of the plastic deformation gradient $\boldsymbol{F}_p$ at the bottom of the melting pool at time $t = 179$ $\mu$s.}
\end{figure}
\begin{figure}[!htb]
\centering
\subfloat[(a)]{\includegraphics[height=0.4\textwidth]{SxxFinal.png}}
\subfloat[(b)]{\includegraphics[height=0.4\textwidth]{SyyFinal.png}}
\newline
\subfloat[(c)]{\includegraphics[height=0.42\textwidth]{SzzFinalM900P900.png}}
\caption{\label{fig:SigmaFpMeltPool} Residual stress components after the cooling stage of the simulation.}
\end{figure}
In order to understand the plastic deformation in the representative volume induced by the laser scan, the components of $\boldsymbol{F}_p$ are shown in Figure \ref{fig:FpMeltPool}(a)-(c) when the laser beam is approximately in the centre of the representative volume.
The thermal expansion ahead of the laser beam leads to plastic compression along the $x$ axis. This is visible on the two sides and at the bottom of the melting pool in Figure \ref{fig:FpMeltPool}(a), where blue areas are present.
By contrast, the plastic deformation along the $y$ axis is compressive at the bottom of the melting pool but it is tensile on the two sides, as shown by the red regions in Figure \ref{fig:FpMeltPool}(b).
Since plastic deformation is isochoric, the region at the bottom of the melting pool expands along the $z$ axis, as shown in Figure \ref{fig:FpMeltPool}(c). Because of the subsequent solidification, the melting pool becomes more shallow and the plastic deformation concentration shown in Figure \ref{fig:FpMeltPool}(a)-(c) spreads in the central region of the representative volume. This is the reason why the regions with the largest plastic deformation are concentrated at depth 30 $\mu$m $< z < 80$ $\mu$m, as shown in Figure \ref{fig:resstress}.
At the top free surface, the plastic deformation is slightly smaller than in the centre of the representative volume, as shown in Figure \ref{fig:schmidfactor}(a)-(c). This is the reason why not all grains with larger Schmid factor at the top surface can accommodate large plastic deformations, as stated in section \ref{sec:grainoripdef}. \\
The diagonal components of the residual Cauchy stress tensor are shown in Figures \ref{fig:SigmaFpMeltPool}(a)-(c). $\boldsymbol{\sigma}_{xx}$ and $\boldsymbol{\sigma}_{yy}$ are heterogeneous at the top surface and in the centre of the representative volume, as shown in Figures \ref{fig:SigmaFpMeltPool}(a) and (b). This shows that they are strongly affected by the grain orientation and that the substrate is applying an important constraint along $x$ and $y$ on the top surface. In fact, even if the top surface is free, large values of $\boldsymbol{\sigma}_{xx}$ and $\boldsymbol{\sigma}_{yy}$ are present. The compressive stress $\boldsymbol{\sigma}_{zz}$ is concentrated in the substrate but it is not very large at the top surface because of the free boundary condition.

\section{Discussion}
\label{sec:discussion}

The simulations suggest the following mechanism for the build up of the residual stresses, as shown in Figure \ref{fig:constraintsubstrate}. The thermal expansion around the molten pool induces plastic compression along the horizontal directions, $x$ and $y$, and plastic tension along the out-of-plane $z$ direction. This is particularly obvious at the bottom of the melting pool, as shown in Figure \ref{fig:FpMeltPool}. After the laser scan is over, tensile residual stresses build up along $x$ and $y$ directions to compensate the plastic compression. This effect is more relevant in the centre of the representative volume, however constraints are present between the regions close to the top surface and the centre. This induces residual stresses along the $x$ and $y$ directions also close to the top surface, as shown in Figure \ref{fig:SigmaFpMeltPool}. \\
\begin{figure}[!htb]
\centering
\includegraphics[width=1.0\textwidth]{ConstraintSubstrate}
\caption{\label{fig:constraintsubstrate} Mechanism for the build up of the residual stresses.}
\end{figure}
The different grain orientations cause the heterogeneity of the residual stress. Grains with slip systems favourably oriented for plastic slip during loads along the $x$, $y$ and $z$ directions are more likely to develop plastic deformation along the corresponding directions. The simulations show that the maximum Schmid factor along the three axes is a suitable quantity to predict the corresponding plastic deformation and consequent residual stress, as shown in Figures \ref{fig:schmidfactor} and \ref{fig:resstress}. \\
The simulations show that, in the stiffness degradation method, the presence of regions with a small residual stiffness imposes a constraint on the regions near the melting pool, preventing more plastic deformation to take place. This is the origin of the underestimation of the plastic deformation and of the residual stress using the stiffness degradation method. By contrast, the element elimination and reactivation method predicts higher plastic deformation in selective laser melting simulations. This difference between the two methods must be taken into account in future studies. \\
Regions with larger tensile residual stresses along specific directions are more likely to undergo fracture perpendicular to those directions. Therefore, the present simulations suggest that fracture perpendicular to the $x$ and $y$ axes is more probable because of the magnitude of the corresponding residual stresses $\boldsymbol{\sigma}_{xx}$ and $\boldsymbol{\sigma}_{yy}$. These two directions have similar residual stresses, as shown in Figure \ref{fig:resstress}. This is consistent with laser-induced microcracking \cite{VRANCKEN2020464}, in which fracture is observed at the top surface perpendicular to the laser scan direction.

\section{Conclusions}
\label{sec:conclusions}

A modelling approach to predict residual stresses at the grain scale in additive manufactured 316 stainless steel is developed.
Thermal-fluid flow and grain growth simulations are used to provide the temperature field and the grain structure respectively, which are implemented in crystal plasticity finite element simulations. A method for element elimination and reactivation is developed to simulate melting and solidification during the laser scan. This method allows to reinitialise state variables, such as the plastic deformation, when elements are reactivated and it represents a step forward compared with previous studies. \\
Simulations show that residual stresses are strongly correlated with the plastic deformation, which in turn depends on the grain orientation. The directions with largest residual stresses are on the laser scan plane, both parallel and perpendicular to the laser scan direction. This is consistent with the microcracks observed experimentally. Even though this study is focused on 316 stainless steel and selective laser melting, the computational method developed can be applied to other metals and other additive manufacturing processes like directed energy deposition and electron beam additive manufacturing. More detailed simulations of residual stresses will be conducted in the future with experimental validation to systematically investigate the relationships between the manufacturing parameters, molten pool flow, grain structures, and residual stresses.

\section*{Acknowledgement}
Thanks to Dr. Yin Zhang from Georgia Institute of Technology for providing EBSD data used to obtain the Euler angles for grain orientations. Thanks to Prof. Yinmin Morris Wang for useful discussion on his experiments on AM 316L stainless steel. This research is supported by the Ministry of Education, Singapore, under its Academic Research Fund Tier 2 (MOE-T2EP50120-0012).

\clearpage

\bibliographystyle{model1-num-names}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}\label{sec:intro}\vspace{-.2cm}
In recent years, there has been a growing trend in data science to develop tools for learning and making inference from signals or data observed on networks. The latter is also known as graph signals which form the subject of investigation in the emerging field of graph signal processing (GSP).
Through modeling real-world networks as graphs and encoding the network data as filtered graph signals, an emerging trend is to develop tools for learning the latent graph topology from these network data; see \cite{dong2019learning, mateos2019connecting} and the references therein.
These tools have widespread applications in the studies of social, financial, and biology networks \cite{newman2018networks}.

Previous works on topology learning with GSP models have focused on mono-layer networks that consist of a single network 'layer', and the tools developed are applicable to single-way data only. However, this may not yield a faithful model for many complex systems since networks and graphs do not simply live in isolation. Instead, many networked systems are better described by multi-layer networks featuring {coupling} interactions \emph{within} and \emph{across} network layers. Examples include: opinion dynamics on correlated topics \cite{parsegov2016novel}, multi-dimensional diffusion \cite{gomez2013diffusion}, protein-protein interactions \cite{prvzulj2011protein}, animal networks \cite{finn2019use}, and relations in image pixels \cite{shi2000normalized}, etc. Observations made on these complex systems are usually generated from two or more coupled networks, and they give rise to multi-attribute observations on nodes, i.e., multi-way graph signals. Naturally, the graph structure embedded in these network data shall be treated using a multilayer graph model as they can not be captured by simple individual networks or flattened networks without structure.

This paper is motivated by the above to develop an inter-layer coupling aware GSP framework for learning an accurate graph structure from multi-way data. As a special case, we focus on modeling a {product graph} \cite{hammack2011handbook} with homogeneous graph layers. We also develop the corresponding inference algorithms.
In light of this, we contribute to both modeling and algorithm aspects in GSP with:
\begin{itemize}[leftmargin=*, noitemsep]
    \item We adopt a general product graph filter for modeling multi-attribute (a.k.a.~multi-way) graph signals. The model uses a bivariate polynomial graph filter to describe interactions over the \emph{physical} and \emph{coupling} graphs. We show that it encompasses a number of common dynamics on multilayer networks.

    \item We develop inference algorithms for topology reconstruction and blind centrality detection in learning the product graph model. Our development relies on an observation that the eigenvectors of graph signals' covariance can be written as Kronecker product of respective eigenvectors of coupling and physical graphs. Leveraging this observation, we suggest spectral methods by proposing two solutions based on the layer/node-wise unfolding and nearest Kronecker product decomposition (NKD). We compare the two solutions by analyzing the effects of inter-layer coupling and show that NKD yields an exact solution under milder assumptions.

\end{itemize}
Lastly, we present numerical experiments on synthetic and real data to corroborate with our analysis. The experiments highlight the effects of inter-layer coupling strengths on topology inference.\vspace{.1cm}

\begin{figure}[t]
\centering
\resizebox{.9\linewidth}{!}{\sf \includegraphics{fig/product_graph/product_final2.png}}\vspace{-.3cm}
\caption{(Left) Physical graph $\mathcal{G}^{\tt G}$ and coupling graph $\mathcal{G}^{\tt C}$, together with the overall product graph ${\cal G}$. (Right) Interaction graph ${\cal G}^{\tt I}$ with different set of parameters $\bm{\gamma}$. }\vspace{-.5cm}
\label{fig:cp}
\end{figure}

\noindent \textbf{Related Works.}
The closest works to ours are \cite{kadambari2021product,kalaitzis2013bigraphical} which studied product graph learning from multi-way graph signals focusing on a smoothness condition defined with the Cartesian product graph. Such conditions restricted the application of these techniques to a special form of coupling mechanism in the multi-way data, which are relaxed by our general product graph filter model.
Besides graph learning, the works \cite{loukas2019stationary,grassi2017time} proposed a time-vertex GSP framework modeling the temporal evolution of graph signals. These models can be regarded as a special case of ours where the coupling graph is fixed as a known path graph.
Moreover, \cite{sandryhaila2014big, stanley2020multiway,zhang2021introducing} studied models for the generic multi-way GSP with possibly heterogeneous graph layers. While the graph model is considerably more generic, these works focused on simple graph filters, which do not account for a number of inter-layer coupling mechanisms. Lastly, we note that \cite{butler2022convolutional} studied an alternative multigraph GSP model without multi-way data.
To our best knowledge, this work introduces the first product graph model with general inter-layer coupling accompanied by a suite of topology inference algorithms.\vspace{-.1cm}

\section{Multi-attribute Graph Signals}\vspace{-.2cm}
This section introduces a generative model for multi-attribute graph signals inspired by network dynamic processes on multiplex graphs with homogeneous graph layers. In particular, our development involves a product graph model \cite{hammack2011handbook} and a general multi-dimensional graph signal/filter model based on the former.\vspace{.1cm}

\noindent \textbf{Product Graph Model.}
The center of our study is a \emph{product graph} ${\cal G} = ({\cal G}^{\tt{C}}, {\cal G}^{\tt{G}})$ formed by two undirected graphs, denoted respectively as $\mathcal{G}^{\tt{C}} =(\mathcal{V}^{\tt{C}}, \mathcal{E}^{\tt{C}},{\bm A}^{\tt{C}})$,  $\mathcal{G}^{\tt{G}}=(\mathcal{V}^{\tt{G}}, \mathcal{E}^{\tt{G}},{\bm A}^{\tt{G}})$ with $|\mathcal{V}^{\tt{C}}|=M$, $|\mathcal{V}^{\tt{G}}|=N$. Note that ${\bm A}^{\tt{C}}$, ${\bm A}^{\tt{G}}$ are weighted adjacency matrices associated with the respective graphs of the edge sets $\mathcal{E}^{\tt{C}} \subseteq \mathcal{V}^{\tt{C}} \times \mathcal{V}^{\tt{C}}$, $\mathcal{E}^{\tt{G}} \subseteq \mathcal{V}^{\tt{G}} \times \mathcal{V}^{\tt{G}}$. We adopt the interpretation of ${\cal G}$ as a \emph{multi-layer graph} such that each node in ${\cal V}^{\tt{G}}$ corresponds to a physical entity, e.g., individual on a social network, etc., and each node in ${\cal V}^{\tt{C}}$ corresponds to a layer/attribute of these physical entities. In this way, ${\cal G}^{\tt{G}}$ will be referred as the \emph{physical graph}, while ${\cal G}^{\tt{C}}$ is the \emph{coupling graph} between layers/attributes; see Fig.~\ref{fig:cp} (left).

Compared to conventional single graph layer models, an important feature of ${\cal G}$ is that it explicitly models \emph{inter-layer coupling} where all $N$ nodes co-exist in the $M$ interacting layers. It is then instrumental to define the \emph{interaction graph} ${\cal G}^{\tt I}$ with the node set ${\cal V}^{\tt I} = {\cal V}^{\tt{C}} \times {\cal V}^{\tt{G}}$. Notice that $|{\cal V}^{\tt I}| = NM$ and the edge set of ${\cal G}^{\tt I}$ can then be encoded into the adjacency matrix:\vspace{-.1cm}
\beq \label{eq:gene_graph}
\A^{\tt I}=\gamma_1\I\otimes \A^{\tt{G}} + \gamma_2\A^{\tt{C}} \otimes\I + \gamma_3\A^{\tt{C}} \otimes\A^{\tt{G}},\vspace{-.1cm}
\eeq
where $\bm{\gamma} = (\gamma_1,\gamma_2,\gamma_3) \geq 0$ satisfying $\sum_{i=1}^3 \gamma_i = 1$ are the \emph{coupling parameters} of the generalized product graph, and $\otimes$ denotes Kronecker product.
Different combinations of $\bm{\gamma}$ can lead to different interaction patterns among the nodes in ${\cal V}$ \cite{tran2018generalized}. For example, with $\bm{\gamma} = (\frac{1}{2},\frac{1}{2},0)$, $\A^{\tt I}$ describes the Cartesian product graph and the matrix is also known as the supra-adjacency matrix; with $\bm{\gamma} = (0,0,1)$, $\A^{\tt I}$ describes the Kronecker product graph; with $\bm{\gamma} = (\frac{1}{3}, \frac{1}{3}, \frac{1}{3})$, $\A^{\tt I}$ describes the strong product graph; see Fig.~\ref{fig:cp} (right).\vspace{.2cm}

\noindent \textbf{Multi-attribute Graph Signals.}
We concentrate on modeling multi-attribute graph signals living on ${\cal G}^{\tt I}$ where the \emph{graph filter} is induced by the graph factors ${\bm A}^{\tt C}, \A^{\tt G}$ and takes an excitation signal as its input. For example, the filter can be a function of ${\bm A}^{\tt I}$.
We consider a general formulation of multi-layer graph filter \cite{natali2020forecasting}:\vspace{-.1cm}
\beq \textstyle \label{eq:gf}
{\cal H}( \A^{\tt C}, \A^{\tt G} ) = \sum_{i=0}^{T_{\tt G}} \sum_{j=0}^{T_{\tt C}} h_{ij} (\A^{\tt{C}})^j \otimes (\A^{\tt{G}})^i,\vspace{-.1cm}
\eeq
where $h_{ij} \in \RR$ are the filter coefficients and $T_{\tt G}, T_{\tt C} \in \ZZ_+ \cup \{ \infty \}$ are the orders of the \emph{bivariate filter polynomial}. The multi-attribute (a.k.a.~multi-way) graph signals are then modeled as:\vspace{-.15cm}
\beq \label{eq:gsp}
{\bm y}^{(s)} = {\cal H}( \A^{\tt C}, \A^{\tt G} ) {\bm x}^{(s)} + {\bm w}^{(s)},\vspace{-.2cm}
\eeq
where $s \in \NN$ denotes the sample index, ${\bm x}^{(s)}, {\bm w}^{(s)} \in \RR^{NM}$ are the excitation signals and observation noise, respectively. Note that ${\bm y}  = ( {\bm y}^{(1)} ; \cdots ; {\bm y}^{(M)} ) \in \RR^{NM}$ whose $m$th block, ${\bm y}^{(m)}$, corresponds to observations on the nodes $\{ (m,i) \in {\cal V}^{\tt I} : i \in {\cal V}^{\tt G} \}$ and the graph signals have been organized in a layer-by-layer fashion.

It is obvious that for all coupling parameters $\bm{\gamma}$, any polynomial of ${\bm A}^{\tt I}$ can be written as \eqref{eq:gf}.
Notably, it is common \cite{gomez2013diffusion} to use the supra-adjacency matrix, ${\bm A}^{\tt I}$ with $\bm{\gamma} = ( \frac{1}{2}, \frac{1}{2}, 0 )$, to describe interactions in a multi-layer graph, and subsequently consider the polynomial of ${\bm A}^{\tt I}$ as the graph filter. However, this is not sufficient to describe certain interactions on ${\cal G}$, e.g., a polynomial of ${\bm A}^{\tt C} \otimes {\bm A}^{\tt G}$ cannot be expressed\footnote{Note that the issue boils down to the non-existence of a polynomial $\tilde{h}(\cdot)$ such that ${h}(ab) = \tilde{h}(a+b)$ for all $h(\cdot), a,b$.} as a polynomial of ${\bm A}^{\tt I}$ with $\bm{\gamma} = ( \frac{1}{2}, \frac{1}{2}, 0 )$.

We conclude this section by showcasing two example dynamic processes to illustrate the generality of our graph filter model:\vspace{-.1cm}
\begin{Exa} Consider the Friedkin-Johnsen multi-dimensional opinion dynamics \cite{friedkin1990social, parsegov2016novel}. At time $t \geq 0$ and the $s$th discussion, the multi-dimensional opinions for  agents in ${\cal V}^{\tt G}$ evolve as \vspace{-.15cm}
\begin{equation} \label{eq:mul_degroot_dyn}
{\bm x}(t+1) = ( {\bm A}^{\tt C} \otimes {\bm A}^{\tt G} ) {\bm x}(t) + {\bm x}^{(s)}\vspace{-.15cm}
\end{equation}
such that the $m$th block of ${\bm x}(t)$ represents the opinions of $N$ agents on the $m$th topic. The weighted adjacency matrices ${\bm A}^{\tt C}, {\bm A}^{\tt G}$ represent the logical dependencies between the $M$ topics, and the mutual trusts between the $N$ agents, respectively. They are properly scaled such that the vector $( {\bm I}_{NM} - {\bm A}^{\tt C} \otimes {\bm A}^{\tt G} ){\bf 1} > {\bm 0}$ represents the self trusts of each agent on the topic. Moreover, ${\bm x}^{(s)}$ is the initial belief of the agents. Under the above premises, we have \vspace{-.15cm}
\beq
{\bm y}^{(s)} = \lim_{t \to \infty} {\bm x}(t) = ({\bm I}_{NM} - {\bm A}^{\tt C} \otimes {\bm A}^{\tt G} )^{-1} {\bm x}^{(s)} , \vspace{-.2cm}
\eeq
where we note that $({\bm I}_{NM} - {\bm A}^{\tt C} \otimes {\bm A}^{\tt G} )^{-1}$ is a special case of \eqref{eq:gf} since the latter is a function of $\A^{\tt I}$ with $\bm{\gamma} = (0,0,1)$.
\end{Exa}
\begin{Exa}
The diffusion process in \cite{gomez2013diffusion} describes dynamics such as social contact, epidemic, etc., on a multi-layer graph. In particular, with the excitation ${\bm x}^{(s)}$, the states of nodes at time $t$ evolve as:\vspace{-.1cm}
\begin{equation} \textstyle
\frac{ d {\bm x}(t) }{dt} = - {\bm x}(t) + \big( \A^{\tt C} \otimes {\bm I}_N + {\bm I}_M \otimes \A^{\tt G} \big) {\bm x}(t) + {\bm x}^{(s)}.\vspace{-.1cm}
\end{equation}
With properly scaled $\A^{\tt C}, \A^{\tt G}$, the above has a unique equilibrium: \vspace{-.15cm}
\begin{equation}
{\bm y}^{(s)} = ( {\bm I}_{NM} - \A^{\tt C} \otimes {\bm I}_N - {\bm I}_M \otimes \A^{\tt G} )^{-1} \x^{(s)}.\vspace{-.1cm}
\end{equation}
where we note that $( {\bm I}_{NM} - \A^{\tt C} \otimes {\bm I}_N - {\bm I}_M \otimes \A^{\tt G} )^{-1}$ is a special case of \eqref{eq:gf} since the latter is a function of $\A^{\tt I}$ with $\bm{\gamma} = (\frac{1}{2},\frac{1}{2},0)$.\vspace{-.1cm}
\end{Exa}
\noindent Graph filter models in the form of \eqref{eq:gf} may also be found in other data scenarios such as the graph causal processes \cite{mei2016signal}. These models may be based on other types of interaction graphs ${\cal G}^{\tt I}$ such as strong product, but they can nevertheless be covered by \eqref{eq:gf}.\vspace{-.2cm}

\section{Product Graph Learning}\vspace{-.2cm}
We aim to infer the product graph model through observing the multi-attribute graph signals in \eqref{eq:gsp}. We concentrate on learning the coupling and physical graphs ${\bm A}^{\tt C}, {\bm A}^{\tt G}$. Our idea is to apply spectral methods with reliable performance regardless of the graph filter \eqref{eq:gf}.

To fix ideas, we denote the eigenvalue decompositions (EVDs) for the adjacency matrices by:\vspace{-.1cm}
\beq
\A^{\tt C} = \V^{\tt C} \bm{\Lambda}^{\tt C} (\V^{\tt C})^\top, \A^{\tt G} = \V^{\tt G} \bm{\Lambda}^{\tt G} (\V^{\tt G})^\top,\vspace{-.1cm}
\eeq
where $\V^{\tt C}, \V^{\tt G}$ are orthogonal and $\bm{\Lambda}^{\tt C}, \bm{\Lambda}^{\tt G}$ are diagonal matrices. Without loss of generality, the eigenvalues are sorted in decreasing order. For example, the $i$th column vector of ${\bm V}^{\tt G}$, denoted ${\bm v}_i^{\tt G}$, corresponds to the $i$th largest eigenvalue ($\lambda_i^{\tt G}$) in ${\bm A}^{\tt G}$. We assume that:\vspace{-.1cm}
\begin{assumption} \label{ass:dis}
For any $i=1,\ldots,N$, $j=1,\ldots,M$, the magnitudes of frequency response $|h(\lambda_j^{\tt C},\lambda_i^{\tt G})|$ have distinct values, where we defined $h(\lambda^{\tt C}, \lambda^{\tt G}) := \sum_{i=0}^{T_{\tt G}}\sum_{j=0}^{T_{\tt C}} h_{ij} (\lambda^{\tt C})^j (\lambda^{\tt G})^i$ according to \eqref{eq:gf}.\vspace{-.1cm}
\end{assumption}
\noindent The above assumption holds for cases when the graph filter is a function of ${\bm A}^{\tt I}$ where ${\bm A}^{\tt C}, {\bm A}^{\tt G}$ have distinct eigenvalues.
Consequently, the EVD of the graph filter \eqref{eq:gf} is derived as (cf.~\cite{sandryhaila2014big}):\vspace{-.1cm}
\beq
{\cal H}( \A^{\tt C}, \A^{\tt G} ) = ( \V^{\tt C} \otimes \V^{\tt G} ) {\cal H}( \bm{\Lambda}^{\tt C}, \bm{\Lambda}^{\tt G} )  ( \V^{\tt C} \otimes \V^{\tt G} )^\top , \vspace{-.1cm}
\eeq
such that ${\cal H}( \bm{\Lambda}^{\tt C}, \bm{\Lambda}^{\tt G} )$ is a diagonal matrix. Under standard white noise conditions for ${\bm x}^{(s)}, {\bm w}^{(s)}$, i.e., both are zero-mean and satisfy $\EE[{\bm x}^{(s)} ({\bm x}^{(s)})^\top] = {\bm I}$, $\EE[{\bm w}^{(s)} ({\bm w}^{(s)})^\top] = \sigma^2 {\bm I}$, the graph signal covariance ${\bm C}_y = \EE[ {\bm y}^{(s)} ({\bm y}^{(s)})^\top ]$ can be derived as:\vspace{-.1cm}
\beq \label{eq:covy}
{\bm C}_y =
( \V^{\tt C} \otimes \V^{\tt G} ) |{\cal H}( \bm{\Lambda}^{\tt C}, \bm{\Lambda}^{\tt G} )|^2  ( \V^{\tt C} \otimes \V^{\tt G} )^\top + \sigma^2 {\bm I}.\vspace{-.1cm}
\eeq
Eq.~\eqref{eq:covy} makes an important observation about the covariance matrix ${\bm C}_y$. When $\sigma = 0$, the eigenvectors of ${\bm C}_y$ are given by the columns of ${\bm V}^{\tt C} \otimes {\bm V}^{\tt G}$. We will demonstrate in \S\ref{sec:estV} that ${\bm V}^{\tt C}, {\bm V}^{\tt G}$ can be retrieved by decomposing the eigenvectors of ${\bm C}_y$.

Suppose for now that
the matrices $\widehat{\bm V}^{\tt C}, \widehat{\bm V}^{\tt G}$ with possibly permuted columns of $\V^{\tt C}, \V^{\tt G}$ are given, we can recover the graph topology and/or detect central nodes in ${\cal G}^{\tt C}, {\cal G}^{\tt G}$ through adapting the spectral methods developed in several existing works. In particular:\vspace{.1cm}

\noindent \textbf{P1. Topology Reconstruction.} It is shown in \cite{segarra2017network} that graph topologies can be reconstructed from the spectral template of graph shift operator (GSO). For example, the following problem recovers ${\bm A}^{\tt G}$:\vspace{-.1cm}
\begin{align}
 \textstyle \min_{ \bm{\lambda}^{\tt G} , \widehat{\bm A}^{\tt G} }~ & \textstyle \|{\rm vec}( \widehat{\bm A}^{\tt G})\|_1 + \frac{\rho}{2} \| \widehat{\bm A}^{\tt G} - \widehat\V^{\tt G} {\rm Diag}(\bm{\lambda}^{\tt G}) (\widehat\V^{\tt G})^\top \|_F^2 \notag \\[.0cm]
\text{s.t.}~ & |{\rm diag}(\widehat{\bm A}^{\tt G})| \leq \epsilon {\bf 1}, \widehat{\bm A}^{\tt G} {\bf 1} \geq {\bf 1}, \label{eq:topo_recon} \\[-.6cm] \notag
\end{align}
where $\rho > 0, \epsilon > 0$ are regularization parameters, and ${\rm vec}(\cdot)$ is the vectorization operator.
Similar formulation can be applied for $\A^{\tt C}$. \vspace{.1cm}

\noindent \textbf{P2. Centrality Estimation.}
It is observed \cite{roddenberry2021blind, he2022detecting, he2021identifying} that eigen-centrality vector can be inferred from graph signals under mild conditions.
Their idea can be adapted for product graphs by noting ${\bm v}_1^{\tt C} \otimes {\bm v}_1^{\tt G}$ is the only positive eigenvector in ${\bm C}_y$. Under the premise that the latter can be decomposed (cf.~\S\ref{sec:estV}), we propose:
\vspace{-.1cm}
\beq \textstyle \label{eq:centrali}
\widehat{\bm c}^{\tt G} = \widehat{\bm v}^{\tt G}_c,~\widehat{\bm c}^{\tt C} = \widehat{\bm v}^{\tt C}_c~~\text{where}~~\widehat{\bm v}_{i^\star} = \widehat{\bm v}^{\tt C}_c \otimes \widehat{\bm v}^{\tt G}_c\vspace{-.1cm}
\eeq
and $i^\star = \argmin_i {\sf P}( \widehat{\bm v}_i )$ with
${\sf P}( {\bm x} ) := \min\{ \| {\bm x} - ({\bm x})^+ \|, \| {\bm x} + (-{\bm x})^+ \| \}$, $({\bm x})^+ := \max\{ {\bf 0}, {\bm x} \}$, and ${\sf P}({\bm x}) = 0$ iff ${\bm x}$ is all-positive or all-negative. Notice $\widehat{\bm v}_i$ denotes the $i$th eigenvector of ${\bm C}_y$.\vspace{-.2cm}

\subsection{Recovering ${\bf V}^{\tt C}, {\bf V}^{\tt G}$ from ${\bf C}_y$} \label{sec:estV} \vspace{-.1cm}
Our remaining task is to recover the eigenvectors in ${\bm V}^{\tt C}, {\bm V}^{\tt G}$ \emph{individually} from observed data ${\bm y}^{(s)}$. For simplicity, we consider the same conditions leading to \eqref{eq:covy} with $\sigma = 0$ and propose two solutions.\vspace{.1cm}

\noindent \textbf{Exact Solution by NKD.} Observe that the diagonal entries in $|{\cal H}(\bm{\Lambda}^{\tt C}, \bm{\Lambda}^{\tt G})|^2$ of \eqref{eq:covy} may not be sorted in descending order. Applying EVD on ${\bm C}_y$ thus produces the eigenvector matrix:
\vspace{-.1cm}
\beq \label{eq:vhat_cov}
\widehat{\bm V} = (\V^{\tt C} \otimes \V^{\tt G}) \bm{\Pi} = \big[ \cdots~{\bm v}_{\pi^{\tt C}(i)}^{\tt C} \otimes {\bm v}_{\pi^{\tt G}(i)}^{\tt G} \cdots \big],\vspace{-.1cm}
\eeq
where $\bm{\Pi}$ is a permutation matrix that orders the columns of $\V^{\tt C} \otimes \V^{\tt G}$ according to the magnitude of frequency response. For any $i \in \{1,\ldots,NM\}$, $\pi^{\tt C}(i), \pi^{\tt G}(i)$ are indices of eigenvectors for $\V^{\tt C}, \V^{\tt G}$ of the $i$th highest frequency response in $|{\cal H}(\bm{\Lambda}^{\tt C}, \bm{\Lambda}^{\tt G})|^2$, i.e., $\widehat{\bm v}_i$.

The structure illustrated in \eqref{eq:vhat_cov} indicates that every column vectors of $\widehat{\V}$ can be written as a Kronecker product.
To this end, we observe the classical result adapted from \cite[Theorem 2.1]{van1993approximation}.\vspace{-.1cm}
\begin{Prop} \label{prop:van}
Consider ${\bm M} \in \RR^{MN \times PQ}$, ${\bm A} \in \RR^{M \times P}$, ${\bm B} \in \RR^{N \times Q}$.
If ${\bm M} = {\bm A} \otimes {\bm B}$, then ${\cal R}({\bm M}) = {\rm vec}( {\bm A} ) {\rm vec}( {\bm B} )^\top$, where \vspace{-.1cm}
\[
{\cal R}( {\bm M} ) = \left[ \begin{array}{c}
{\bm M}_1 \\ \vdots \\ {\bm M}_P
\end{array} \right],~
{\bm M}_i = \left[
\begin{array}{c}
{\rm vec}({\bm M}_{1,i})^\top \\ \vdots \\ {\rm vec}({\bm M}_{M,i})^\top
\end{array} \right],~i=1,\ldots,P.\vspace{-.1cm}
\]
Note ${\bm M}$ has been partitioned into $N \times Q$ blocks and ${\bm M}_{i,j}$ is the $(i,j)$th block.
Subsequently, ${\bm A}, {\bm B}$ can be uniquely  recovered (up to a scalar factor) through a suitable decomposition of ${\cal R}({\bm M})$.\vspace{-.1cm}
\end{Prop}
\noindent When ${\bm M} \neq {\bm A} \otimes {\bm B}$, e.g., with noisy observations, \cite{van1993approximation} consider the \emph{nearest Kronecker product decomposition} (NKD) problem:\vspace{-.1cm}
\beq \textstyle \label{eq:nkp}
\begin{array}{rl}
\min_{ {\bm A}, {\bm B} , \alpha } & \| {\cal R}( {\bm M} ) - \alpha \, {\rm vec}( {\bm A} ) {\rm vec}( {\bm B} )^\top \|_F^2 \\
\text{s.t.} & \| {\rm vec}( {\bm A} ) \| = 1,~\| {\rm vec}( {\bm B} ) \| = 1,\vspace{-.1cm}
\end{array}
\eeq
which can be solved by finding the top singular vectors of ${\cal R}( {\bm M} )$.

Proposition~\ref{prop:van} shows that ${\bm v}_{\pi^{\tt C}(i)}^{\tt C}, {\bm v}_{\pi^{\tt G}(i)}^{\tt G}$ can be recovered by applying NKD on ${\cal R}( \widehat{\bm v}_i )$ from \eqref{eq:vhat_cov}. Consequently, collecting the NKD outputs on ${\cal R}( \widehat{\bm v}_i )$ for all $i$ forms two matrices whose columns are permuted and repeated copies of the columns of $\V^{\tt C}, \V^{\tt G}$. This suggests applying the Gram-Schmidt procedure on these matrices [cf.~\eqref{eq:gs}] for estimating $\V^{\tt C}, \V^{\tt G}$; see Algorithm~\ref{algo:prod_learn}.
The algorithm returns an exact solution under H\ref{ass:dis} and other conditions:
\vspace{-.1cm}
\begin{Corollary}
Under H\ref{ass:dis} and noiseless observations ($\sigma=0$) with $S \to \infty$, Algorithm~\ref{algo:prod_learn} recovers the columns of $\V^{\tt C}, \V^{\tt G}$.\vspace{-.1cm}
\end{Corollary}
\noindent We remark that for {\bf P2}, one may skip step 4 of Algorithm~\ref{algo:prod_learn} to obtain $i^\star = \argmin_i {\sf Pos}( \widehat{\bm v}_i^{\sf noisy} )$ in \eqref{eq:centrali} and thus the pair $( \widehat{\bm v}_{i^\star}^{\tt C}, \widehat{\bm v}_{i^\star}^{\tt G} )$. In this way, H\ref{ass:dis} can be further weakened to guarantee exact recovery.

\begin{figure}[t]
\vspace{-.3cm}
\begin{algorithm}[H]
\caption{Learning Product Graph ({\bf P1} and/or {\bf P2})}
\label{algo:prod_learn}
\begin{algorithmic}[1]
\STATE \textbf{INPUT}: Set of $S$ observed graph signals $\{ {\bm y}^{(s)} \}_{s=1}^S$.
\STATE Evaluate $\CovYhat^S = (1/S)\sum_{s=1}^S {\bm y}^{(s)} ({\bm y}^{(s)})^\top$ and compute its eigenvectors as $\hatV^{\sf noisy} = {\sf EVD}( \CovYhat^S )$.
\STATE For any $i=1,\ldots,NM$, solve the NKD problem for ${\cal R}( \widehat{\bm v}_i^{\sf noisy} )$ [cf.~\eqref{eq:nkp}] to obtain the pair $( \widehat{\bm v}_i^{\tt C}, \widehat{\bm v}_i^{\tt G} )$.
\STATE Perform Gram-Schmidt (GS) to obtain the orthogonal matrices:\vspace{-.1cm}
\beq \label{eq:gs}
\widehat{\V}^{\tt C} = {\sf GS}( [\widehat{\bm v}_1^{\tt C}, \ldots, \widehat{\bm v}_{NM}^{\tt C}]),~
\widehat{\V}^{\tt G} = {\sf GS}([\widehat{\bm v}_1^{\tt G}, \ldots, \widehat{\bm v}_{NM}^{\tt G}]).\vspace{-.45cm}
\eeq
\STATE Apply the methods in {\bf P1} for topology reconstruction \eqref{eq:topo_recon} or {\bf P2} for centrality detection \eqref{eq:centrali}.
\end{algorithmic}
\end{algorithm}
\vspace{-3em}
\end{figure}

\vspace{.15cm}
\noindent \textbf{Simplified Solution by Unfolding.} We conclude the section by proposing and analyzing an alternative to NKD in Algorithm~\ref{algo:prod_learn} for estimating $\V^{\tt C}, \V^{\tt G}$.
The alternative solution is inspired by \cite{sandryhaila2014big,zhang2021graph} through unfolding the multi-attribute graph signals into layer-wise and node-wise matrices. Interestingly, we show that this simplified design can be as effective as NKD, but only when the frequency response of the product graph filter \eqref{eq:gf} is \emph{separable} \cite{grassi2017time}.

\renewcommand{\UrlFont}{\normalsize}

To this end, we denote the unfolding of ${\bm y}^{(s)}$ in \eqref{eq:gsp} as ${\bm Y}^{(s)} = [ {\bm y}_1^{(s)}, \ldots, {\bm y}_M^{(s)} ]$ and consider ${\bm C}_y^{\sf layer} = \EE[ ({\bm Y}^{(s)})^\top {\bm Y}^{(s)} ]$, ${\bm C}_y^{\sf node} = \EE[ {\bm Y}^{(s)} ({\bm Y}^{(s)})^\top ]$. Their resultant eigenvector matrices $\widetilde{\bm V}^{\tt C} = {\sf EVD}( {\bm C}_y^{\sf layer} )$, $\widetilde{\bm V}^{\tt G} = {\sf EVD}( {\bm C}_y^{\sf node} )$ can then be used in lieu of the estimates in step~4 of Algorithm~\ref{algo:prod_learn}. The following proposition analyzes the covariance matrices:\vspace{-.15cm}
\begin{Prop} \label{prop:sep}
Under H\ref{ass:dis}. Assume that the excitation ${\bm x}^{(s)}$ is zero-mean satisfying $\EE[ {\bm x}^{(s)} ( {\bm x}^{(s)} )^\top ] = {\bm I}$ and ${\bm w}^{(s)} = {\bm 0}$. It holds:\vspace{-.2cm}
\beq
\begin{split}
{\bm C}_y^{\sf layer} & \textstyle = \sum_{i=1}^M {\bm v}_i^{\tt C} ( {\bm v}_i^{\tt C} )^\top \sum_{j=1}^N | h ( \lambda_i^{\tt C}, \lambda_j^{\tt G} )|^2, \\
{\bm C}_y^{\sf node} & \textstyle = \sum_{j=1}^N {\bm v}_j^{\tt G} ( {\bm v}_j^{\tt G} )^\top \sum_{i=1}^M | h ( \lambda_i^{\tt C}, \lambda_j^{\tt G} )|^2.\\[-.15cm]
\end{split}
\eeq
\end{Prop}
\noindent The proof can be found in the appendix.
To gain insight, we focus on node-wise covariance ${\bm C}_y^{\sf node}$ whose eigenvectors $\widetilde{\V}^{\tt G}$ correspond to desired $\V^{\tt G}$ iff the eigenvalues $\sum_{i=1}^M | h ( \lambda_i^{\tt C}, \lambda_j^{\tt G} )|^2$ are distinct. However, this scenario of repeated eigenvalues is common for ${\bm C}_y^{\sf node}$ when $h ( \lambda_i^{\tt C}, \lambda_j^{\tt G} ) = h ( \lambda_i^{\tt C} \lambda_j^{\tt G} )$, i.e., the graph filter is based on interaction graphs with $\bm{\gamma} = (0,0,1)$. For example, when $h ( \lambda_i^{\tt C}, \lambda_j^{\tt G} ) = e^{ \lambda_i^{\tt C} \lambda_j^{\tt G} }$ with the spectrum for ${\cal G}^{\tt C}$ as $\{-1,1,2\}$, then ${\bm C}_y^{\sf node}$ has repeated eigenvalues if $\{ 0, -1.27 \}$ are in the spectrum of ${\cal G}^{\tt G}$.

On the other hand, a sufficient condition for ${\bm C}_y^{\sf node}$ to admit distinct eigenvalues hinges on the \emph{separable filter} property \cite{grassi2017time}.
Concretely, we need $h(\lambda^{\tt C}, \lambda^{\tt G}) = h^{\tt C}(\lambda^{\tt C}) h^{\tt G}( \lambda^{\tt G})$ for all $\lambda^{\tt C}, \lambda^{\tt G}$. For example, when $h ( \lambda_i^{\tt C}, \lambda_j^{\tt G} ) = e^{ \frac{1}{2}(\lambda_i^{\tt C} + \lambda_j^{\tt G}) }$, the interaction graph is free of the direct Kronecker product. In this case, the eigenvectors of ${\bm C}_y^{\sf layer}, {\bm C}_y^{\sf node}$ serve as good surrogates for $\V^{\tt C}, \V^{\tt G}$, respectively.\vspace{-.1cm}

\vspace{-.2cm}
\section{Numerical Experiments}\vspace{-.2cm}

\noindent \textbf{Synthetic Data.}
We consider graph filters that are based on ${\bm A}^{\tt I}$ in \eqref{eq:gene_graph} with the parameters $\bm{\gamma} = ( \gamma_1, 2 \gamma_1, 1-3\gamma_1 )$. As $\gamma_1 \downarrow 0$, the interaction graph has strong \emph{inter-layer coupling} of Kronecker product form.
We fix $\mathcal{G}^{\tt C}$ to be the tree graph in Fig.~\ref{fig:cp} with $M=3$ where ${\bm A}^{\tt C}$ is the unweighted adjacency matrix; while ${\bm A}^{\tt G}$ will be the unweighted adjacency matrix for ${\cal G}^{\tt G}$ to be determined below.\vspace{.1cm}

\noindent {\textbf{Example 1: Topology Reconstruction.}}
This example examines the performance of product graph topology reconstruction ({\bf P1}). We generate ${\cal G}^{\tt G}$ as an Erdos-Renyi graph with connection probability of $p=0.4$. We benchmark the proposed `NKD' based Algorithm~\ref{algo:prod_learn} against `Unfold' which estimates $\widetilde{\V}^{\tt C}, \widetilde{\V}^{\tt G}$ from the layer-wise/node-wise unfolded graph signals; `PGL' refers to the product graph learning method based on smoothness of layer-wise/node-wise unfolded graph signals in \cite{kadambari2021product}; `Flatten' refers to applying \cite{segarra2017network} directly to  ${\bm C}_y$. For the {\tt SpecTemp} problem \eqref{eq:topo_recon} in Algorithm~\ref{algo:prod_learn}, `Unfold', and `Flatten', we take $\rho = 40$, $\epsilon = 10^{-6}$ and solve \eqref{eq:topo_recon} by \texttt{cvx}.

We consider synthetic data generated from \eqref{eq:gsp} with the graph filter given by $ \HAJ=e^{\tau \A^{\tt I}}$, where $\A^{\tt I}$ is parameterized by $\gamma_1 \in [0,\frac{1}{3}]$ as described. We also set $\tau = 1 / \max_i d_i$ with $d_i = \sum_{j=1}^N \A^{\tt I}_{ij}$. Excitation and noise signals satisfy ${\bm x}^{(s)} \sim {\cal N}( {\bm 0}, {\bm I} )$, ${\bm w}^{(s)} \sim {\cal N}( {\bm 0}, 0.01 {\bm I} )$.
Fig.~\ref{fig:topology} compares the performance of topology reconstruction in terms of the $F_1$ score against the sample size $S$ and graph size $N$. We assume known $\gamma_1$ and evaluate the $F_1$ scores by comparing the ground truth to $\widehat{\A}^{\tt I}$ reconstructed by \eqref{eq:gene_graph} using the estimated $\A^{\tt C}, \A^{\tt G}$.
We observe that `NKD' recovers the graph topology for a wide range of sample sizes and graph sizes, regardless of the coupling parameter $\gamma_1$. Meanwhile, `Unfold' and `PGL' are sensitive to $\gamma_1$ -- their performances are comparable to Algorithm~\ref{algo:prod_learn} when $\gamma_1 = \frac{1}{3}$, but it degrades noticeably when $\gamma_1 = 0.01$. In all cases, `Flatten' fails to estimate the graph topology.\vspace{.1cm}

\begin{figure}[t]
\centering
\resizebox{0.995\linewidth}{!}{\sf \input{fig/simulation/topology/new_simu/topo_withoutD_new}}\vspace{-.4cm}
\caption{\textbf{Topology Reconstruction.} $F_{1}$ score against (Left) sample size $S$ with $N =10$; (Right) graph size $N$ with $S=1000$.
}\vspace{-.4cm}
\label{fig:topology}
\end{figure}

\noindent \textbf{Example 2: Central Nodes Detection.}
The next example examines the performance of central node detection ({\bf P2}). We focus on detecting the central nodes of ${\cal G}^{\tt G}$ using \eqref{eq:centrali}, where ${\cal G}^{\tt G}$ is generated as a core-periphery graph \cite{core_per}. The node set $\mathcal{V}^{\tt G}=\{1,\ldots,N\}$ is partitioned into $\Vcore=\{1,\ldots,10\}$, $\Vpe=\mathcal{V}^{\tt G}\backslash\Vcore$. For every pair $(i,j) \in \mathcal{V}^{\tt G} \times \mathcal{V}^{\tt G}$, edges are assigned randomly with probability $1$ if $i,j \in \Vcore$; with probability $0.2$ if $i \in \Vcore, j \in \Vpe$; with probability $0.05$ if $i,j \in \Vpe$. Our goal is to detect the $10$ nodes in $\Vcore$ based on the observed multi-attribute graph signals. We focus on comparing the proposed `NKD' based Algorithm~\ref{algo:prod_learn} with `Unfold' which estimates $\widetilde{\V}^{\tt C}, \widetilde{\V}^{\tt G}$ from the node-wise unfolded graph signals.

We consider synthetic data modeled after \eqref{eq:gsp} with two graph filters: (a) ${\cal H}^{\sf inv}( \A^{\tt C}, \A^{\tt G} ) =(\I- \tau_1 \A^{\tt I})^{-1}$, (b) ${\cal H}^{\sf exp}( \A^{\tt C}, \A^{\tt G} ) = e^{\tau_2 \A^{\tt I}}$ with $\tau_1 = 1/\max_i d_i$, $\tau_2 = 10\tau_1$. We set ${\bm x}^{(s)} \sim {\cal N}( {\bm 0}, {\bm I} )$, ${\bm w}^{(s)} \sim {\cal N}( {\bm 0}, 0.01 {\bm I} )$.
Fig.~\ref{fig:cen_de} compares the error rate in identifying the central nodes $\Vcore$ as top-10 central nodes \eqref{eq:centrali}. As seen, `NKD' delivers lower error rate with smaller number of samples for the strong coupling scenario ($\gamma_1 = 0.01$). The result corroborates with Proposition~\ref{prop:sep} as the covariance with node/layer-wise unfolded signals have close (or repeated) eigenvalues, significantly increasing the sample complexity in eigenvector estimation.\vspace{.1cm}

\begin{figure}[t]
\centering
\resizebox{0.995\linewidth}{!}{\sf \input{fig/simulation/centrality/cen_de}}\vspace{-.4cm}
\caption{\textbf{Central Nodes Detection.} Error rate against (Left) sample size $S$ with $N=80$; (Right) graph size $N$ with $S=5MN$.
}\vspace{-.2cm}
\label{fig:cen_de}
\end{figure}

\begin{figure}[t]
\centering
{\resizebox{0.47\linewidth}{!}{\sf \includegraphics{fig/simulation/rollcall/GraphG_revise.png}}
}
{\resizebox{0.45\linewidth}{!}{\sf \includegraphics{fig/simulation/rollcall/GC_revise.png}}
}%
\vspace{-.25cm}
\caption{\textbf{US Senate Rollcalls Data.} (Left) Estimated Senate topology ${\cal G}^{\tt G}$ color coded by the two communities found by spectral clustering on ${\bm A}^{\tt G}$.
(Right) Estimated coupling graph of the topics of rollcalls.
}\vspace{-.55cm}
\label{fig:rollcall}
\end{figure}

\noindent \textbf{Real Data: US Senate Roll Calls.}
We apply Algorithm~\ref{algo:prod_learn} to reconstruct the topology ({\bf P1}) from the 113th-116th US Senate rollcalls [available: \url{https://voteview.com/data}].
We first apply NMF \cite{arora2012learning} to perform topic modeling
on rollcalls' descriptions to get $M=4$ topics, we then model and infer a product graph with $M=4$ layers and $N=50$ nodes, where each node corresponds to the Senators from a state.
Fig.~\ref{fig:rollcall} shows the learnt graph topology from $S=101$ samples through applying Algorithm~\ref{algo:prod_learn}. The estimated topology is reasonable: the state graph ${\cal G}^{\tt G}$ identifies the Republican and Democratic states as two clusters, and the coupling graph ${\cal G}^{\tt C}$ shows a strong connection between topics on `circuit judges' and `fiscal'.\vspace{.1cm}

\noindent \textbf{Conclusions.} This paper considered the problem of learning from multi-attribute graph signals by proposing a general product graph filter model and developing its inference algorithms.
Future works include analyzing the sampling complexity of Algorithm~\ref{algo:prod_learn} and studying the effects of layer-wise/node-wise unfolding of graph signals on sampling complexity or identifiability.

\newpage

\bibliographystyle{IEEEtran}