\section{Introduction}

In porous metals, void coalescence often drives the onset of the macroscopic flow localisation that marks the end of uniform deformation and acts as a precursor to failure, as well as the initiation and propagation of ductile cracks \cite{tekouglu2015localization, guo2018void, liu2019micromechanism}. Previous studies suggest that for conventional plasticity theory, where no material length scale enters the constitutive law (absence of stress/strain gradient induced size effect), a decrease in the inter-void spacing promotes void coalescence \cite{pardoen2000extended, srivastava2013void} and results in the collapse of the yield surface \cite{torki2015void, torki2017theoretical}. While for a fixed inter-void spacing, it is well established that the imposed stress state has a pronounced effect on the onset of void coalescence in conventional plasticity theory. For example, it has been shown that an increase in the imposed stress triaxiality (a ratio of the first to second stress invariant) promotes void growth and early onset of void coalescence \cite{budiansky1982void, koplik1988void, needleman1992void, benzerga2010ductile}. Void coalescence is simply the event where the plastic flow localises within the inter-void ligaments and successively links the neighboring voids \cite{koplik1988void}. The plastic flow localisation within the inter-void ligament, however, will induce plastic strain gradients that in turn may affect the strengthening and hardening of the material. This raises a fundamental question: how does the interaction of inter-void spacing (or ligament size) and the gradient induced material size effect, affect the localisation of plastic flow causing void coalescence for a given stress state?

The gradient induced size effect resulting in strengthening and hardening in metals has been confirmed in many material tests involving non-uniform deformation including indentation \cite{stelmashenko1993,ma1995size}, torsion \cite{fleck1994}, and bending \cite{stolken1998}. The size dependent material response on the micron scale in metal plasticity implies that the growth of micron sized voids also exhibits significant size effects \cite{tvergaardniordson_2007, niordson_2008}. In general, it has been shown that the gradient induced size effect leads to slower growth rates for smaller voids \cite{tvergaard2004nonlocal, li2006rve, monchiet2013gurson, nielsen_2016, holte2019}. An accurate representation of void coalescence due to plastic flow localisation within micron sized inter-void ligaments, therefore, also requires material models that represent stresses over the relevant length scales. Phenomenological theories describing the strengthening and hardening due to plastic strain gradients express the plastic work in terms of both plastic strain and plastic strain gradient, thereby introducing a length scale into the material model. Herein, the strain gradient plasticity theory proposed by Gudmundson \cite{gudmundson2004} is used, which includes both dissipative (non-recoverable) and energetic (recoverable) gradient contributions within a small strain formulation based on visco-plasticity. The mathematical formulation and associated variational structure originate from Fleck and Willis \cite{fleck2009}, and the material model is implemented into the commercial finite element software ABAQUS using a user element (UEL) subroutine \cite{martinez2019}.

The objective of this work is to understand the interaction of the inter-void spacing (or ligament size) and the resultant gradient induced material size effect on void coalescence for a range of imposed stress states. To achieve this, three dimensional finite element unit cell calculations for a periodic array of initially spherical voids embedded in a strain gradient enhanced material matrix are carried out. Several unit cell geometries have been analyzed to investigate the effect of inter-void ligament size under multiple loading conditions. The imposed stress states are characterised by fixed values of the stress triaxiality and the Lode parameter (a measure of the third stress invariant). The value of the Lode parameter is shown to affect the evolution of voids in computations involving conventional plasticity theory \cite{zhang2001numerical, kim2004modeling, gao2006modeling, barsoum2007rupturemicro, srivastava2013void} and in experiments \cite{bao2004fracture, barsoum2007rupture, srivastava2012effect} only at relatively low stress triaxiality levels. However, it is likely that in an anisotropic material matrix \cite{srivastava2015effect} with anisotropy introduced by the void distribution \cite{srivastava2013void}, as for the present investigation, the effect of the Lode parameter can be important even at high stress triaxialities.

Our results show that for a conventional material matrix, increasing the inter-void ligament size results in an increase in the critical stress to void coalescence, up to a threshold value of inter-void ligament size. The sensitivity of the critical stress to the inter-void ligament size is found to increase with increasing stress triaxiality. The quantitative effect of the Lode parameter is found to be small for the stress triaxiality values varying from $1$ to $3$. However, for inter-void ligament sizes below the threshold value, the critical stress is smallest for a Lode parameter value of $-1$, whereas above the threshold value the critical stress is smallest for a Lode parameter value of $0$. For a void in a strain gradient enhanced material matrix, the value of the critical stress for void coalescence increases with increasing length parameter i.e. increasing gradient effect. This effect of the length parameter on the critical stress magnitude is found to increase with increasing imposed stress triaxiality and decreasing inter-void ligament size. This is because at higher stress triaxiality values and for smaller inter-void ligament sizes, there is an increase in the propensity for plastic flow localisation that introduces strong plastic strain gradients and in turn hardens the ligament. This mechanism leads to a decrease in the dependence of critical stress on the inter-void ligament size with increasing length parameter. The gradient induced strengthening also tends to homogenize the deformation in the unit cell thus decreasing the effect of the Lode parameter.

The structure of the manuscript is as follows. Section~\ref{S:problemformulation} frames the study and presents the numerical method. The unit cell geometries considered, the method utilized to impose proportional loading throughout the deformation history, and the strain gradient plasticity material model are also presented in Section \ref{S:problemformulation}. The numerical results are presented and discussed in Section \ref{S:results}. Finally, the key results and conclusions of this work are summarized in Section \ref{sec:Conclusion}.

\section{Problem formulation and modelling approach}
\label{S:problemformulation}

This work considers a limit load-type analysis to determine the critical stress level at which a given microstructure configuration loses load carrying capacity. Hence, an elastic-perfectly plastic material model is employed. The configuration of the unit cell and the simulation setup is described in Section \ref{sec:geometryandmesh}, while the approach to prescribe a constant value of stress triaxiality and Lode parameter is outlined in Section \ref{sec:numericalmethod}.

\subsection{Unit cell geometry and FE mesh}
\label{sec:geometryandmesh}

Three dimensional finite element calculations are carried out to model the response of an array of spherical voids with initial radius $r_0$, Fig.~\ref{fig:cell}. The unit cell has edge lengths $2a_i^0$ along the three coordinate axes, $x_i$ ($i=1,2,3$), and inter-void spacings thereby equal $2l_i^0=2a_i^0-2r_0$. Symmetry about three planes perpendicular to the coordinate axes implies that only $1/8$ of the unit cell needs to be modelled.

\begin{figure}[htbp!]
\centerline{\includegraphics[width=3.5in]{Cell-setup.png}}
\caption{Schematic showing the periodic arrangement of voids in the $x_2$- and $x_3$-plane. The distribution along the $x_1$-direction is not shown for simplicity.}
\label{fig:cell}
\end{figure}

For all unit cells considered, the initial void volume fraction is $f_0=0.01$, where $f_0=(4/3\pi r_0^3)/(8a^0_1 a^0_2 a^0_3)$. The initial void radius, $r_0$, is kept constant, while the cell dimensions are varied to achieve various initial inter-void spacings as in ref.~\cite{srivastava2013void}. The geometric parameters for the different cell dimensions are given in Table~\ref{tab:geometry}. For all unit cells, $a^0_1/r_0=a^0_2/r_0$. Finite element meshes for four unit cell configurations are shown in Fig.~\ref{fig:meshes}. The modelling setup does not account for softening due to void evolution since a small strain formulation is used. It is assumed that the unit cell represents the material condition immediately before failure, neglecting the deformation history leading to this state. Hence, model predictions for the loss of load carrying capacity signal the onset of localisation. The critical equivalent stress at the onset of localisation is recorded and reported in the results section.

\begin{table}[htbp!]
\centering
\caption{Geometric parameters for the various unit cells considered for $f_0=0.01$. Based on \cite{srivastava2013void}.}
\label{tab:geometry}
\begin{tabular}{cccccc}
\hline
$a^0_1/r_0=a^0_2/r_0$ & $a^0_3/r_0$ & $l^0_1/r_0=l^0_2/r_0$ & $l^0_3/r_0$ \\
\hline
$6.06$   & $1.43$        & $5.05$        & $0.43$        \\
$5.55$   & $1.70$        & $4.55$        & $0.70$        \\
$5.21$   & $1.94$        & $4.21$        & $0.95$        \\
$4.97$   & $2.12$        & $3.97$        & $1.12$        \\
$4.58$   & $2.50$        & $3.58$        & $1.50$        \\
$4.18$   & $3.00$        & $3.18$        & $2.00$        \\
$3.75$   & $3.75$        & $2.75$        & $2.75$        \\
\hline
\end{tabular}
\end{table}

\begin{figure}[htb!]
\begin{overpic}[width=\linewidth]{Meshes-shade.png}
\put(18,51){(a)}
\put(70,51){(b)}
\put(18,0){(c)}
\put(70,0){(d)}
\end{overpic}
\caption{Finite element meshes showing 1/8 of the unit cell with an initially spherical void of radius $r_0$ in the centre giving an initial void volume fraction of $f_0=0.01$ for (a) $\nicefrac{l_1^0}{r_0}=\nicefrac{l_2^0}{r_0}=5.06$; $\nicefrac{l_3^0}{r_0}=0.43$, (b) $\nicefrac{l_1^0}{r_0}=\nicefrac{l_2^0}{r_0}=4.21$; $\nicefrac{l_3^0}{r_0}=0.95$, (c) $\nicefrac{l_1^0}{r_0}=\nicefrac{l_2^0}{r_0}=3.58$; $\nicefrac{l_3^0}{r_0}=1.5$ and (d) $\nicefrac{l_1^0}{r_0}=\nicefrac{l_2^0}{r_0}=2.75$; $\nicefrac{l_3^0}{r_0}=2.75$. The number of elements ranges from 1896 for (a) to 2512 for (d).}
\label{fig:meshes}
\end{figure}

\subsection{Numerical method}
\label{sec:numericalmethod}

The unit cells are subject to prescribed displacements and the boundary conditions applied to the faces of the cell are

\begin{align}
\label{eq:BC}
    u_1(a_1^0, x_2, x_3)=U_1(t),\quad T_2(a_1^0, x_2, x_3)=T_3(a_1^0, x_2, x_3)=0 \nonumber \\
    u_2(x_1, a_2^0, x_3)=U_2(t),\quad T_1(x_1, a_2^0, x_3)=T_3(x_1, a_2^0, x_3)=0 \nonumber \\
    u_3(x_1, x_2, a_3^0)=U_3(t),\quad T_1(x_1, x_2, a_3^0)=T_2(x_1, x_2, a_3^0)=0
\end{align}

\noindent The applied symmetry boundary conditions are

\begin{align}
\label{eq:BC_symm}
    u_1(0, x_2, x_3)=0, \quad T_2(0, x_2, x_3)=T_3(0, x_2, x_3)=0 \nonumber \\
    u_2(x_1, 0, x_3)=0, \quad T_1(x_1, 0, x_3)=T_3(x_1, 0, x_3)=0 \nonumber \\
    u_3(x_1, x_2, 0)=0, \quad T_1(x_1, x_2, 0)=T_2(x_1, x_2, 0)=0
\end{align}

In Eq. \eqref{eq:BC}, $U_1(t)$ is prescribed and the time history of the displacements $U_2(t)$ and $U_3(t)$ are determined such that a prescribed stress state is maintained. The loading direction is fixed in stress space by enforcing constant ratios between the normal stress components throughout the deformation history such that

\begin{align}
\label{eq:rho}
    \Sigma_{22} = \rho_2\Sigma_{11}, \qquad \Sigma_{33} = \rho_3\Sigma_{11},
\end{align}

\noindent where $\rho_2$ and $\rho_3$ are constants. The overall stress components $\Sigma_{ij}$ are found by volume averaging over all elements, such that: $\Sigma_{ij}=\int_V\sigma_{ij}\text{d}V/V$, where $V$ is the unit cell volume.

The overall effective stress, $\Sigma_e$, and the overall hydrostatic stress, $\Sigma_h$, are given by

\begin{align}
    \Sigma_e=\frac{1}{\sqrt{2}}\sqrt{(\Sigma_{11}-\Sigma_{22})^2+(\Sigma_{22}-\Sigma_{33})^2+(\Sigma_{33}-\Sigma_{11})^2}, \nonumber\\ \qquad \Sigma_h=\frac{1}{3}(\Sigma_{11}+\Sigma_{22}+\Sigma_{33})\nonumber,
\end{align}
which in terms of the relative stress ratios become
\begin{align}
    \Sigma_e=\Sigma_{11}\frac{1}{\sqrt{2}}\sqrt{(1-\rho_2)^2+(\rho_2-\rho_3)^2+(\rho_3-1)^2}, \\ \Sigma_h=\Sigma_{11}\frac{1}{3}(1+\rho_2+\rho_3).
\end{align}
The stress triaxiality, $T$, and the Lode parameter, $L$, are given by
\begin{align}
    T=\frac{\Sigma_h}{\Sigma_e}=\frac{\sqrt{2}}{3}\frac{1+\rho_2+\rho_3}{\sqrt{(1-\rho_2)^2+(\rho_2-\rho_3)^2+(\rho_3-1)^2}}
    \label{eq:T}
\end{align}
and
\begin{align}
    L=\frac{2\Sigma_{22}-\Sigma_{11}-\Sigma_{33}}{\Sigma_{11}-\Sigma_{33}}=\frac{2\rho_2-1-\rho_3}{1-\rho_3}.
    \label{eq:L}
\end{align}
The overall effective strain, $\overline{E}_e$, is given by
\begin{align}
    \overline{E}_e=\frac{\sqrt{2}}{3} \sqrt{(E_{11}-E_{22})^2+(E_{22}-E_{33})^2+(E_{33}-E_{11})^2}
\end{align}
where the strain components, $E_{ij}$, are found in a way analogous to the stress components.

\subsubsection{Multiple point constraints}

The macroscopic normal stress components vary throughout the deformation, but the stress ratios are maintained in each increment of the simulation according to Eq.~\eqref{eq:rho}. This is achieved by creating multi point constraints through the user subroutine MPC in ABAQUS, which enables enforcing relationships between degrees of freedom in one or more nodes.

Additional degrees of freedom are added to impose boundary conditions on all sides of the model while prescribing the stress ratios. Three dummy nodes, $N_i$, are created outside of the mesh and connected to one connector node, $M$, in the mesh as shown in Fig.~\ref{fig:MPC}. This connection is made through spring elements (SPRING2 elements from the ABAQUS element library). The displacement in the $x_1$-direction is then prescribed at the $N_1$-dummy node, while the displacements (in the $x_2$- and $x_3$-directions) corresponding to the desired stress triaxiality and Lode parameter are calculated and applied to the $N_2$- and $N_3$-dummy nodes. The displacement of the connector node, $M$, is coupled to the displacement of the nodes located at $(a_1^0,x_2,x_3)$, $(x_1,a_2^0,x_3)$ and $(x_1,x_2,a_3^0)$ in the direction of the respective face normals. In this way, the displacement of the dummy nodes, $N_i$, is linked to the unit cell.

\begin{figure}[htb!]
\centering
\centerline{\includegraphics[width=2.5in]{MPC.PNG}}
\caption{The spring elements for the multiple point constraints connected to one connector node, $M$, in the finite element mesh.}
\label{fig:MPC}
\end{figure}

The displacement of the dummy nodes, $N_i$, is related to the forces, $F_i$, at the faces of the unit cell through

\begin{align}
   F_i=k_i(u_i^{N_i}-u_i^M) \quad \text{with} \quad i=1,2,3,
\label{eq:springs}
\end{align}

\noindent and $k_i$ being the spring element constants given by $k_i=E(\nicefrac{A_i}{a_i^0})\times10^{-1}$, where the factor of $10^{-1}$ is introduced to stabilise the numerical solution, following Ref.~\cite{tekoglu2014representative}. The forces, $F_i$, are the resultant of all traction across the corresponding surface and relates to the macroscopic stresses through

\begin{align}
    \Sigma_{11}=\frac{F_1}{A_1}, \qquad A_1=a_2^0a_3^0 \nonumber \\
    \Sigma_{22}=\frac{F_2}{A_2}, \qquad A_2=a_1^0a_3^0 \nonumber \\
    \Sigma_{33}=\frac{F_3}{A_3}, \qquad A_3=a_1^0a_2^0
\label{eq:stress}
\end{align}

\noindent where $A_i$ is the area over which the forces act. Combining Eqs.~\eqref{eq:rho}, \eqref{eq:springs}, and \eqref{eq:stress}, gives the dummy node displacements

\begin{align}
    \rho_2=\frac{\Sigma_{22}}{\Sigma_{11}}=\text{const.} \Rightarrow u_2^{N_2}=u_2^M+\rho_2\frac{A_2}{A_1}\frac{k_1}{k_2}(u_1^{N_1}-u_1^M) \nonumber \\
    \rho_3=\frac{\Sigma_{33}}{\Sigma_{11}}=\text{const.} \Rightarrow u_3^{N_3}=u_3^M+\rho_3\frac{A_3}{A_1}\frac{k_1}{k_3}(u_1^{N_1}-u_1^M),
\end{align}

\noindent where $\rho_2$ and $\rho_3$ are input values for the stress ratio, $u_i^{N_j}$ is the displacement of dummy node $j$ in the direction of $x_i$, $u_i^M$ is the displacement in $x_i$-direction of the connector node $M$, $A_i$ are areas from Eq.~\eqref{eq:stress}, and $k_i$ are the spring element constants. Another relevant procedure for imposing multiple point constraints without spring elements can be found in \cite{liu2016void}.

The calculations are carried out for three values of Lode parameter, $L=-1, 0$, and $1$. The Lode parameter values $L=-1$ ($\Sigma_{11}>\Sigma_{22}=\Sigma_{33}$) and $L=1$ ($\Sigma_{11}=\Sigma_{22}>\Sigma_{33}$) correspond to overall axisymmetric stress states, while $L=0$ ($\Sigma_{11}>\Sigma_{22}>\Sigma_{33}$) correspond to an overall state of shear plus hydrostatic stress. For each value of Lode parameter, three triaxialities are considered, $T=1, 2$, and $3$. The values for $\rho_2$ and $\rho_3$ to achieve these stress states are given in Table~\ref{tab:input}.

\begin{table}[htb!]
\caption{Input parameters determining the prescribed stress state.}
\label{tab:input}
\begin{tabular}{m{0.1\textwidth}m{0.1\textwidth}m{0.1\textwidth}m{0.1\textwidth}}
\hline
L & T & $\rho_2$ & $\rho_3$  \\
\hline
\multirow{4}{*}{-1} & $\frac{1+2\rho_2}{3(1-\rho_2)}$ & $\rho_3$ & $\frac{1+2\rho_2}{3(1-\rho_2)}$ \\
                    & 1 & 0.4 & 0.4    \\
                    & 2 & 0.625 & 0.625  \\
                    & 3 & 0.727273 & 0.727273  \\
\hline
\multirow{4}{*}{0} & $\frac{1+\rho_2}{\sqrt{3}(1-\rho_2)}$ & $\frac{1+\rho_3}{2}$ & $\frac{\sqrt{3}T-1}{\sqrt{3}T+1}$ \\
                   & 1 & 0.634 & 0.268 \\
                   & 2 & 0.776 & 0.552 \\
                   & 3 & 0.8386 & 0.6772 \\
\hline
\multirow{4}{*}{1} & $\frac{2+\rho_2}{3(1-\rho_2)}$  &  1 &  $\frac{3T-2}{3T+1}$ \\
                   & 1 & 1 & 0.25 \\
                   & 2 & 1 & 0.57 \\
                   & 3 & 1 & 0.70 \\
\hline
\end{tabular}
\end{table}

The calculations were carried out using the commercial finite element code ABAQUS with the gradient theory applied to the matrix material through a UEL subroutine. The reader is referred to Ref. \cite{martinez2019} for further details on the implementation. The calculations use 20-node user defined elements. The number of elements in the finite element meshes is varied from a minimum of 1896 to a maximum of 2512 elements, Fig.~\ref{fig:meshes}.

\subsection{Material model: Strain gradient plasticity}
\label{sec:gradientmodel}

The gradient enhanced constitutive model employed is based on the visco-plastic strain gradient plasticity theory proposed by Gudmundson \cite{gudmundson2004} in the context of the mathematical formulation in terms of minimum principles proposed by Fleck and Willis \cite{fleck2009}. For the dissipative version considered, the theory accounts for internal elastic energy storage due to elastic strain and dissipation due to the plastic strain rate, $\dot{\varepsilon}^p_{ij}$, and its spatial gradient, $\dot{\varepsilon}^p_{ij,k}$. Contributions from plastic strain gradients to free energy is ignored. The Principle of Virtual Work (PVW) in Cartesian components is expressed by

\begin{multline}
    \int_{V}\left(\sigma_{ij}\delta\dot{\varepsilon}_{ij}+(q_{ij}-s_{ij})\delta\dot{\varepsilon}^p_{ij}+\tau_{ijk}\delta\dot{\varepsilon}^p_{ij,k}\right)\text{d}V= \\ \int_{S}\left(T_i\delta\dot{u}_i+t_{ij}\delta\dot{\varepsilon}^p_{ij}\right)\text{d}S
\end{multline}

\noindent where $\sigma_{ij}$ and $s_{ij}=\sigma_{ij}-\frac{1}{3}\delta_{ij}\sigma_{kk}$ are the Cauchy stress tensor and the stress deviator, respectively. The micro-stress, $q_{ij}$, is work conjugate to the plastic strain rate, $\dot{\varepsilon}^p_{ij}$, and $\tau_{ijk}$ is a higher order stress, work conjugate to the plastic strain rate gradient, $\dot{\varepsilon}^p_{ij,k}$. The right hand side of the PVW includes the conventional traction, $T_{i}=\sigma_{ij}n_j$ work conjugate to the boundary displacement rate, $\dot u_i$, and the higher order traction, $t_{ij}=\tau_{ijk}n_k$, work conjugate to the plastic strain rate, $\dot{\varepsilon}_{ij}^p$. Here, the outward unit normal to the surface $S$ is $n_i$.  Balance laws for the stress quantities are given by

\begin{gather}
    \sigma_{ij,j}=0 \quad \text{and} \quad q_{ij}-s_{ij}-\tau_{ijk,k}=0
\end{gather}

\noindent where, the first set of equations is the conventional equilibrium equations in the absence of body forces, and the second set is the higher order equilibrium equations. The higher order boundary conditions are imposed such that the void surface is higher order traction free, while symmetry conditions are imposed at the exterior of the cell through $\varepsilon_{12}=0$.

\subsection{Constitutive equations}

The rate-dependent visco-plastic formulation employs a potential to account for plastic dissipation as follows

\begin{equation}
  \Phi\left[\dot{E}^p,E^p\right]=\int_{0}^{\dot{E}^p}\sigma_c\left[\dot{E}^{p'},E^p\right]\text{d}\dot{E}^{p'}
\end{equation}

Here, $\sigma_c$ is the gradient enhanced effective stress, related to the current matrix flow stress through $\sigma_c=\sigma_F[E^p]\left(\frac{\dot{E}^p}{\dot{\varepsilon}_0}\right)^m$, with $\dot{\varepsilon}_0$ denoting the reference strain rate, and $m$ denoting the rate-sensitivity exponent. The material in this work does not undergo strain hardening, making $\sigma_F$ independent of $\dot{E}^p$ and equal to the material yield stress $\Sigma_0$. The viscoplastic law is implemented following the algorithm presented in Ref. \cite{FuentesVisco} to efficiently approach the rate-independent limit. A gradient enhanced effective plastic strain rate is given by

\begin{equation}
\label{eq:dotE}
    \left(\dot{E}^p\right)^2 = \frac{2}{3}\dot{\varepsilon}^p_{ij}\dot{\varepsilon}^p_{ij}+L_D^2\dot{\varepsilon}^p_{ij,k}\dot{\varepsilon}^p_{ij,k}
\end{equation}

\noindent and the associated work conjugate gradient enhanced effective stress by

\begin{equation}
\label{eq:sigmac}
    \sigma^2_c=\frac{3}{2}q_{ij}q_{ij}+\frac{1}{L_D^2}\tau_{ijk}\tau_{ijk}.
\end{equation}

\noindent Here, $L_D$ is a dissipative constitutive length parameter that enters for dimensional consistency. The superscript $D$ refers to dissipative quantities, and the dissipative stress quantities are given by

\begin{equation}
\label{eq:disstress}
    q^D_{ij}=\frac{2}{3}\sigma_c\frac{\dot{\varepsilon}^p_{ij}}{\dot{E}^p}, \,\,\,\,\,\,\,\, \tau^D_{ijk}=L_D^2\sigma_c\frac{\dot{\varepsilon}^p_{ij,k}}{\dot{E}^p}.
\end{equation}

The dissipative length parameter controls the strengthening size effect with an increase in the dissipative length parameter giving an increase in the apparent yield stress in the presence of strain gradients, see~\cite{martinez2016finite,voyiadjis2012thermo}. This work is a limit load analysis, which, by definition, is done to determine the overall yield criterion for a given, specific configuration. Limit load analyses normally idealise materials as perfectly plastic. To avoid strain hardening from the energetic gradient contributions, the energetic length parameter, $L_E$, has been set to zero in this work, and, consequently, the corresponding energetic quantities are omitted.

\section{Numerical results and discussion}
\label{S:results}

Throughout, the following material parameters are used; $\Sigma_{0}/E=0.001$, $\nu=0.3$ and $m=0.01$, where $\Sigma_{0}$ is the yield stress, $E$ is Young's modulus, $\nu$ is the Poisson ratio, and $m$ is the strain rate sensitivity exponent. The value of $m$ is considered sufficiently small for the results to approximate a rate-independent material response. The influence of the Lode parameter, $L$, the stress triaxiality, $T$, and the normalised length parameter, $L_D/r_0$, is studied. The effect of the inter-void ligament size is discussed in combination with the other parameters, $L$, $T$, and $L_D/r_0$.

\subsection{Critical equivalent stress at localisation}
\label{sec:criticalstress}

Figure~\ref{fig:stress-strain} presents the equivalent stress-strain curves for two distinct Lode parameters, $L=-1$ and $1$, but for a fixed stress triaxiality, $T=3$, and a fixed inter-void ligament size of $l_3/r_0=1.5$. The equivalent stress-strain curves are depicted for three length parameters, being, $L_D/r_0=0.2, 0.5$, and $1$ as well as for the conventional limit where $L_D/r_0=0$.

The material response shows a clear effect of plastic strain gradients, such that the larger the length parameter, the higher the equivalent stress level. This means that an increase in the stress level is obtained when down-scaling the microstructure and, thus, yielding of the material is delayed due to increasing strain gradient strengthening. The critical equivalent stress, $\Sigma_e^c/\Sigma_0$, signaling localisation (and coalescence) is taken to be at the plateau of the equivalent stress-strain curve. Several ways exists to establish a coalescence criteria based on either critical stress or strain. The method employed in this work is inspired by the work of~\cite{tekoglu2012criterion}. The critical stresses have been extracted from the end of the equivalent stress-strain curves (as shown in Figures 4a and b). Taking the example of Fig.~\ref{fig:stress-strain}a, the conventional material ($L_D/r_0=0$), the difference between the equivalent stress at $E_e=0.02$ and $0.01$ is less than $0.04\%$.

\begin{figure*}[htb!]
\begin{minipage}[c]{0.45\textwidth}
    \begin{overpic}[width=\textwidth]{Stress-strainL-1.png}
    \put(17,12){\textbf{$L=-1$}}
    \end{overpic}
\subcaption{}
\end{minipage}
\hspace{0.2cm}
\begin{minipage}[c]{0.45\textwidth}
    \begin{overpic}[width=\textwidth]{Stress-strainL1.png}
    \put(17,12){\textbf{$L=1$}}
    \end{overpic}
\subcaption{}
\end{minipage}
\caption{Equivalent stress-strain curve for an inter-void ligament size of $l_3/r_0=1.5$ under loading conditions giving Lode parameters of (a) $L=-1$ and (b) $L=1$ and a triaxiality of $T=3$.}
\label{fig:stress-strain}
\end{figure*}

\subsection{Conventional material: Effect of the inter-void ligament size}
\label{s:Lodeparameter}

The conventional limit, $L_D/r_0=0$, is considered to set the scene for the study of material size effects. The focus here is the effect of inter-void ligament size on the critical stress at localisation under various loading conditions.

First, three values of the Lode parameter are considered, $L=-1, 0$, and $1$, for a fixed stress triaxiality, $T=2$. Figure~\ref{fig:Lode} shows the critical equivalent stress as a function of the inter-void ligament size. For the six smallest inter-void ligaments, the critical equivalent stress is seen to increase when the inter-void ligament becomes bigger irrespective of the value of the Lode parameter. The increase in the critical stress ties to localisation occurring more easily in small inter-void ligaments lowering the load carrying capacity of the unit cell. As the inter-void ligament size increases, the $l_3$-ligament can sustain a higher stress level before localisation, leading to an increase in critical equivalent stress. Also, for the six smallest inter-void ligaments, an increase in the critical stress is found with increasing Lode parameter values. Thus, the lowest critical equivalent stress is found for $L=-1$. The dependence on the Lode parameter can be rationalised by considering the imposed stress state. In comparison to the other cases, the relative stress component, $\rho_3$, is the largest when $L=-1$ (see Tab.~\ref{tab:input}), and localisation is therefore expected in the $l_3$-ligament at a lower overall deformation. In contrast, the $\rho_3$ takes the lowest value for $L=1$, resulting in delayed localisation and the highest critical equivalent stress obtained. In Ref. \cite{srivastava2013void}, void coalescence was found to occur along the ligament with the smallest applied stress for $L>-1$. For all Lode parameters, the relative stress component in the $l_3$-ligament will be smallest as $\rho_3$ is always the lowest stress ratio. For $L=-1$, coalescence occurs in the direction of the smallest inter-void ligament size. This corresponds to the $l_3$-ligament for all geometries except when $l_3/r_0=2.75$ as this is a perfect cube, Table~\ref{tab:geometry}.

\begin{figure}[htb!]
\centerline{\includegraphics[width=3.34in]{Lode.png}}
\caption{Critical equivalent stress vs. normalized inter-void ligament size for three values of the Lode parameter with $T=2$ and $L_D/r_0=0$.}
\label{fig:Lode}
\end{figure}

There is, however, a shift in the localisation pattern when the $l_3$-ligament becomes sufficiently wide, for example, a drop in the coalescence stress is found for $l_3/r_0=2.75$ for $L=0$. The load carrying capacity of the material increases when the distribution of voids diverge from a regular array, i.e. when $l_3/r_0 \neq 2.75$. However, the load carrying capacity will decrease if the arrangement of the voids is such that the inter-void ligament size, in any direction, is too small (e.g. $l_3/r_0 <  2$). At the configuration with $l_3/r_0=2.75$, there is no bias towards the $l_3$-ligament since the unit cell takes a cubic shape. The shift in the localisation is especially prominent for $L=0$ (a state of combined hydrostatic tension and shear) where plastic flow localises at $\approx 45^\circ$ across the cubic unit cell leading to an early loss of load carrying capacity. The shift in the localisation is demonstrated by depicting the contours of the effective plastic strain for two distinct unit cells ($l_3/r_0=1.5$ and $2.75$) subjected to $L=0$ and $T=2$ in Figs.~\ref{fig:contour_Case5+7} and \ref{fig:Shearband}. The material response remains conventional such that $L_D/r_0=0$, and the loading conditions are described by $L=0$ and $T=2$. For the conventional material, the second term of Eq.~\eqref{eq:dotE} is zero ($L_D=0$) and the term \textit{gradient enhanced effective plastic strain} refers to the time integration of only the first term of Eq. \eqref{eq:dotE}. For the elongated unit cell ($l_3/r_0=1.5$), localisation is seen to occur in the smallest ligament, $l_3$, whereas localisation is seen to occur along two corners of the cubic unit cell ($l_3/r_0=2.75$), indicating that the deformation is localised along $\approx 45^\circ$ i.e. across the diagonal. Figure~\ref{fig:contour_Case5+7} shows the contour of equivalent plastic strain across the faces of the cell, while Fig.~\ref{fig:Shearband}(a) and (b) show the contour of the effective plastic strain in the diagonal cross-section of both unit cells at an overall effective strain of $\overline{E}_e=0.03$. By comparing the two contours it is seen that plastic flow is observed across the entire cross-section indicating localisation at $\approx 45^\circ$ for the cubic model, $l_3/r_0=2.75$. In contrast, the plastic flow is constricted for $l_3/r_0=1.5$ shown top right in Fig.~\ref{fig:Shearband}.

\begin{figure}[htb!]
\centering
\begin{overpic}[width=\linewidth]{Case5+7_L0T2.PNG}
\end{overpic}
\caption{Distribution of effective plastic strain for $L=0$, $T=2$, $L_D/r_0=0$, for $l_3/r_0=1.5$ to the left and $l_3/r_0=2.75$ to the right at a macroscopic effective strain of $\overline{E}_e=0.03$. For $l_3/r_0=1.5$, localisation is favoured in the smallest ligament, $l_3$. For the cubic unit cell, however, there is no bias towards any of the ligaments and deformation localises along $\approx 45^\circ$, i.e. across the diagonal.}
\label{fig:contour_Case5+7}
\end{figure}

\begin{figure}[htbp!]
\centering
\begin{overpic}[width=0.5\textwidth]{Case5+7_shear.png}
\put(0,35){(a)} \put(0,0){(b)}
\end{overpic}
\caption{Distribution of effective plastic strain along a cut from corner to corner for $L=0$, $T=2$, $L_D/r_0=0$ at $\overline{E}_e=0.03$ for two geometries: (a) $l_3/r_0=1.5$ and (b) $l_3/r_0=2.75$. In (a) plasticity has not localized along $45^\circ$ and for this geometry localisation is favoured in the smaller ligament, while (b) shows that a band, indicated by the dotted line, has formed at a $45^o$ angle to the main loading axis (the $x_1$-axis), thus lowering the critical effective stress.}
\label{fig:Shearband}
\end{figure}

Next, the effect of stress triaxiality is considered. In Fig.~ \ref{fig:triaxiality}, the critical stress for a conventional material, $L_D/r_0=0$, is shown as a function of the inter-void ligament size for $T=1, 2,$ and $3$ for a fixed value of the Lode parameter, $L=-1$. In the conventional limit, a high level of stress triaxiality yields low critical stress for all ligament sizes considered. The reason being that a high stress triaxiality corresponds to higher relative stress components, $\rho_2$ and $\rho_3$. Figure~\ref{fig:triaxiality} shows little effect of ligament size on the critical equivalent stress for the low value of triaxiality. For $T=1$, the relative stress transverse to the main loading direction is insufficient to invoke localisation in the inter-void ligament and the effect of the ligament size will be limited. The cell instead undergoes macroscopic localisation and, consequently, does not exhibit a profound dependence on the inter-void ligament size. This is in line with results presented in Ref. \cite{tekouglu2015localization}, where $T=1$ has been found to be the limit below which the onset of macroscopic localisation is essentially simultaneous with void coalescence. The results for $T=2$ and $T=3$ in Fig.~\ref{fig:triaxiality} show that the critical equivalent stress is dependent on inter-void ligament size.

A small drop in the critical equivalent stress is seen to occur for the largest ligament for all values of triaxiality in Fig.~\ref{fig:triaxiality}. The effect is most prominent for the highest triaxiality, $T=3$. Figure \ref{fig:triaxiality_symmetry} shows the contour gradient enhanced effective plastic strain of the cubic cell ($l_3/r_0=2.75$) at an effective stress of $\overline{E}_e=0.08$. At a sufficiently large strain, plasticity is seen to initiate at the corner opposite to the void. Due to the symmetry of both the loading condition ($\Sigma_2=\Sigma_3$ for $L=-1$) and the unit cell, bands of plastic deformation are observed to stretch across the $x_1-x_2$ and $x_1-x_3$ faces, ultimately lowering the coalescence stress giving the drop as seen in Fig.~\ref{fig:triaxiality}.

\begin{figure}[htb!]
\centering
    \includegraphics[width=3.34in]{Triaxiality.png}
    \caption{Critical equivalent stress vs. normalized inter-void ligament size for three values of the stress triaxialities with $L=-1$ and $L_D/r_0=0$.}
    \label{fig:triaxiality}
\end{figure}

\begin{figure}[htb!]
\centering
    \includegraphics[width=2.75in]{Triaxiality_symmetry.png}
    \caption{Bands of plastic flow in the cubic unit cell ($l_3/r_0=2.75$) and $L_D/r_0=0$ at an overall equivalent strain of $\overline{E}_e=0.08$. The loading conditions applied to give an axisymmetric stress state with $L=-1$ and $T=2$. Note the rotated coordinate system to show the symmetry of the plastic flow given by the cubic unit cell and $\rho_2=\rho_3$ for $L=-1$.}
    \label{fig:triaxiality_symmetry}
\end{figure}

\subsection{Gradient enriched material: Effect of the inter-void ligament size}
\label{s:lengthscaleparamter}

The effect of gradient strengthening in the matrix material is introduced through the length parameter $L_D$ (see Section~\ref{sec:gradientmodel}). One can imagine down-scaling the microstructure when increasing the value of $L_D/r_0$. Three values of the length parameter, $L_D/r_0=0.2, 0.5$, and $1$, are considered in the following for all combinations of the Lode parameter, $L=-1, 0, 1$, and stress triaxiality, $T=1, 2, 3$. Figure~\ref{fig:SigC_l3} shows the critical effective stress, $\Sigma_e^c/\Sigma_0$, as a function of the inter-void ligament size, $l_3/r_0$, for all combinations. The results obtained for a conventional material, $L_D/r_0=0$ are presented as a reference (see Section~\ref{s:Lodeparameter}).

\begin{figure*}
\centering
\begin{overpic}[width=\textwidth]{Matrix.PNG}
\put(6.2,56.5){(a)} \put(37,56.5){(b)} \put(67,56.5){(c)}
\put(6.2,33){(d)} \put(37,33){(e)} \put(67,33){(f)}
\put(6.2,10){(g)} \put(37,10){(h)} \put(67,10){(i)}
\end{overpic}
\caption{The critical equivalent stress as a function of the smallest inter-void ligament size. Three values of the Lode parameter are considered, $L=-1, 0$, and $1$. For each Lode parameter, three values of the stress triaxiality are considered, $T=1, 2$, and $3$. Throughout, the parameters $\Sigma_0/E=0.001$, $\nu=0.3$ and $m=0.01$ are used. The initial void volume fraction is, $f_0=0.01$. The length parameter that enters through the gradient plasticity theory is $L_D/r_0=0.2, 0.5$ and $1$. A conventional material is modelled with $L_D/r_0=0$ and used as a reference.}
\label{fig:SigC_l3}
\end{figure*}

The general observation is that the critical stress at localisation increases with the magnitude of the length parameter (down-scaling the microstructure). However, the critical stress has a natural upper bound where the gradient strengthening is so severe that the entire matrix material yields. At such large values of the length parameter, the effect of the Lode parameter, stress triaxiality, and inter-void ligament size vanish and the critical equivalent stress is identical for all combinations of geometry and loading condition. The threshold value is evident from Figs.~\ref{fig:SigC_l3}(a)-(i).

Figure~\ref{fig:contour_lengthscale}(a)-(c) show how the length parameter affects the plastic flow in the unit cells by comparing contours of the gradient enhanced plastic equivalent strain, $E^p$, (see Eq.~\eqref{eq:dotE}) for a ligament size of $l_3/r_0=1.5$ subject to $L=-1$ and $T=3$. The contours are extracted at an overall equivalent strain of $\overline{E}_e=0.02$. Figure~\ref{fig:contour_lengthscale}(a) displays the conventional material response where localisation occurs in the $l_3$-ligament. At the same level of the overall deformation, a significantly lower effective plastic strain is observed in Fig.~\ref{fig:contour_lengthscale}(b) and (c) when increasing the length parameter. For $L_D/r_0=0.2$, Fig.~\ref{fig:contour_lengthscale}(b), some plasticity is seen to develop in the $l_3$-ligament, but far less than in the conventional case, while the plasticity has barely initiated at this level of the deformation for $L_D/r_0=0.5$ Fig.~\ref{fig:contour_lengthscale}(c). The corresponding equivalent stress is shown in Figs.~\ref{fig:contour_lengthscale}(d)-(f). It is seen that the level of stress in the unit cell increases with increasing length parameter. The critical equivalent stress is seen to increase with increased length parameter in Fig.~\ref{fig:SigC_l3} and the material can therefore withstand higher stresses.

\begin{figure}
\centering
\begin{overpic}[width=0.5\textwidth]{Contours_stress_strain.png}
\put(18,67){(a)} \put(55,67){(b)} \put(82,67){(c)}
\put(18,-2){(d)} \put(55,-2){(e)} \put(82,-2){(f)}
\end{overpic}
\vspace*{-3mm}
\caption{Distribution of gradient enhanced effective plastic strain for $L=-1$, $T=3$, $l_3/r_0=1.5$ for (a) the conventional material, $L_D/r_0=0$, (b) $L_D/r_0=0.2$ and (c) $L_D/r_0=0.5$ at a macroscopic effective strain of $\overline{E}_e=0.02$. The effective stress, $\Sigma_e$, for the same configuration is shown in the cells at the bottom, also here with (d) the conventional material, (e) $L_D/r_0=0.2$ and (f) $L_D/r_0=0.5$.}
\label{fig:contour_lengthscale}
\end{figure}

Figure~\ref{fig:deformationmode} shows the change in the deformation mechanism that occurs with increased gradient strengthening. The contour of the normalised rate of equivalent plastic strain, $\dot{E}^p/\dot{\overline{E}}_e$, is shown for the unit cell with the smallest inter-void ligament, $l_3/r_0=0.43$ under loading conditions giving $L=1$ and $T=3$ for the conventional material with $L_D/r_0=0$ and the material with the greatest gradient strengthening contribution, $L_D/r_0=1$. Figure~\ref{fig:deformationmode}(a) shows that plastic deformation has developed and localised in the $l_3$-ligament, as expected for a conventional material at this loading condition. For the matrix surrounding the $l_3$-ligament, plasticity is reduced in favour of localisation in the $l_3$-ligament. However, for the gradient strengthened material, plasticity is not only less developed, in line with the gradient strengthening, but also smeared out across the unit cell, see Fig.~\ref{fig:deformationmode}(b). Localisation is to a little extent observed in the $l_3$-ligament, but overall the entire cell experiences plasticity. This is indicative of a change in deformation mechanism along the lines of the one observed in Ref. \cite{tekouglu2015localization} for a stress triaxiality of 1. However, here it is seen with an increasing length parameter. As $L_D/r_0$ increases, the cell is more likely to undergo simultaneous macroscopic localisation and void coalescence in contrast to a conventional material where the cell predominantly undergoes void coalescence for the same loading conditions.

\begin{figure}[htbp!]
\centering
\begin{overpic}[width=0.5\textwidth]{Deformationmodes.png}
\put(2,11){(a)} \put(54,11){(b)}
\end{overpic}
\caption{Change in deformation mode with increased length parameter for $l_3/r_0=0.43$ with loading conditions described by $L=1$ and $T=3$. The conventional material with $L_D/r_0=0$ is shown in a), while b) shows a gradient enriched material with $L_D/r_0=1$.}
\label{fig:deformationmode}
\end{figure}

The combined effect of the stress triaxiality (for a fixed Lode parameter) and the length parameter is visualised by the rows in Fig.~\ref{fig:SigC_l3}, while the combined effect of the Lode parameter (for a fixed stress triaxiality) and length parameter is visualised by the columns in Fig.~\ref{fig:SigC_l3}. Qualitatively, for a fixed stress triaxiality value, the length parameter has a nearly identical impact for all values of the Lode parameter; the critical stress increases with increasing length parameter. It is, however, interesting that the drop in coalescence stress the cubic unit cell ($l_3/r_0=2.75$) subject to $L=0$ diminishes with increasing length parameter for all values of stress triaxiality, Figs.~\ref{fig:SigC_l3}(d)-(f). This is because increased gradient strengthening delays the intensification of the plastic flow and homogenizes the plastic strain field.

For the lowest stress triaxiality value, $T=1$, the effect of the length parameter is small. Nonetheless, the plastic strain gradients that build up around the void give rise to the small increase in gradient strengthening. For the loading conditions giving $T=1$, the onset of localisation is significantly delayed, thus allowing the material to withstand higher critical stress with a smaller dependence on the inter-void ligament size.  The deformation mechanism prevailing at this low value of triaxiality, where macroscopic localisation and void coalescence occur simultaneously \cite{tekouglu2015localization}, implies that the gradients surrounding the void will not influence the critical equivalent stress to a great extent, as the deformation takes place in the entire unit cell.

For a higher value of stress triaxiality, $T=2$, the effect of the length parameter is more prominent as seen in Fig.~\ref{fig:SigC_l3}, and the smaller the inter-void ligament size, the greater the effect of the length parameter is. This is because, at higher stress triaxiality values, the plastic flow tends to localise in the inter-void ligaments as the ligaments diminish in size. The localisation induces large plastic strain gradients that in turn contribute to strengthening. The gradient induced strengthening in the inter-void ligament then inhibits further plastic flow localisation and delays void coalescence. Although not shown here, for $L=-1$ and $T=2$, the gradient strengthening is sufficiently large that increasing the value of $L_D/r_0$ from $1$ to $2$, has a negligible effect. For the intermediate length parameter, $L_D/r_0=0.5$, the effects of triaxiality and inter-void ligament size are still visible but greatly reduced due to the smaller degree of gradient strengthening. For the smallest value of the length parameter, $L_D/r_0=0.2$, the critical equivalent stress values follow those of the conventional material, just at a higher relative level for all inter-void ligament size considered. For $L=0$ and $L=1$, for $T=2$, the same effect is seen. The most pronounced effect of the length parameter is seen for $L=-1$, $T=3$ and $l_3/r_0=0.43$ as this configuration has the lowest critical effective stress for the conventional material, but shows the same critical stress for $L_D/r_0=1$ as in the remaining results.

\section{Summary and conclusions}
\label{sec:Conclusion}
The interaction of the inter-void ligament size and the gradient induced material size effect on void coalescence is investigated for a range of imposed stress states, here characterised by fixed values of the stress triaxiality and the Lode parameter. To this end, three dimensional finite element unit cell calculations for a single initially spherical void embedded in strain gradient enhanced material matrix are carried out. A conventional material matrix (absence of gradient induced strengthening effects) is considered as reference. Increasing the length parameter, and thereby the gradient effect, is equivalent to down-scaling the microstructure. All microstructures considered in this work contain voids that are below a critical flaw size \cite{martinez2019role}. Thus, plasticity theory will reign the material response as the voids are considered too small to be treated as cracks.

The results for the conventional material show that the critical coalescence stress increases when increasing the inter-void ligament size. The effect of the inter-void ligament size is, however, dependent on the imposed stress triaxiality, such that the effect of the inter-void ligament size increases with increasing stress triaxiality. However, above a certain threshold for the inter-void ligament size, the results show a slight decrease in the critical stress. This drop has to do with a transition from plastic flow localisation within the smallest inter-void ligament to plastic flow localization at $\approx 45^\circ$ to the main loading axis. The transition in the plastic flow localisation pattern is found to be particularly pronounced for a Lode parameter of $L=0$. However, irrespective of the Lode parameter value, the transition occurs as the unit cells approach a cubic geometry.

For a void embedded in a strain gradient enhanced material matrix, the value of the critical coalescence stress increases with increasing length parameter i.e. increasing the gradient strengthening effect. The effect of the length parameter is found to intensify with increasing imposed stress triaxiality and decreasing inter-void ligament size. This is due to a propensity for plastic flow localisation in the inter-void ligament when the ligament is small and the stress triaxiality high. Plastic flow localisation introduces large plastic strain gradients which in turn strengthens the ligament and delays further localisation of plastic flow. The strengthening from plastic strain gradients also leads to a weakened dependency in the critical coalescence stress on the inter-void ligament size.  Finally, the results show that there exists a natural upper bound where the gradient strengthening is so severe that the entire matrix material yields. For very large values of the length parameter, the effect of the imposed stress state and the inter-void ligament size vanish, and the critical equivalent stress is identical for all combinations of the unit cell geometry and the loading conditions considered.

\begin{acknowledgment}
\noindent This research was financially supported by the Danish Council for Independent Research through the research project ``Advanced Damage Models with InTrinsic Size Effects'' (Grant no: DFF-7017-00121). The computing resources provided by the high performance research computing center at Texas A\&M University are gratefully acknowledged.
\end{acknowledgment}

\bibliographystyle{asmems4}


\section{Theoretical considerations}

	In the theoretical part of the Supplementary Material, we give the proofs for the propositions in the main text. We first introduce the geometric representation of a qudit state. Exploiting this representation, we prove the disk conjecture (Theorem 1) regarding the geometric form of the maximal maskable set for $d$-dimensional qudit states, and subsequently prove the inverse proposition in qubit case, that a two-dimensional quantum disk must also be maskable. We also explicitly derive the form of isometry for qubit masking. Finally, we prove (Theorem 2) that any qudit broadcastable set is a proper subset of some qudit maskable sets. Moreover, the qudit maskable set has nonzero measure in $d$-dimensional Euclidean space.

	\subsection{Geometric representation of qudit}

	An arbitrary qudit state $\rho$ can be expressed as follows using the $SU(d)$ basis $\{\Lambda_i\}^{d^2-1}_{i=1}$, where $\Lambda_i$ are orthogonal eigenstates satisfying $\Lambda_i^+=\Lambda_i$, $\mathrm{Tr}(\Lambda_i)=0$, and $\mathrm{Tr}(\Lambda_i\Lambda_j)=2\delta_{ij}$. Namely\,\cite{Mahler95, Kimura03},

	\begin{eqnarray}
		\rho=\frac{1}{d}I_d+\frac{1}{2}\sum^{d^2-1}_{i=1}x_i \Lambda_i,~~~ \sum^{d^2-1}_{i=1}x^2_i\leq r_d^2,~~~ x_i\in \mathbb{R},
	\end{eqnarray}
	where $r_d$ satisfies $${\frac{2}{d(d-1)}}\leq r_d^2 \leq {\frac{2(d-1)}{d}}.$$

	The definition of the hyperdisk is crucial for the comprehension of quantum information masking. A hyperdisk $\mathcal{D}$ is defined as the intersection of a hypersphere and a hyperplane in the Euclidean space. Namely,

	\begin{align}
		\mathcal{D}=\{(x_1,x_2,...,x_d):\sum^d_{i=1}x^2_i\leq r^2_d,\sum^d_{i=1}a_ix_i+b=0,a_i \mathrm{~not~ all~ be~ zero}\}.
	\end{align}

	\subsection{General form of maskable set in arbitrary dimensions}

	\textit{Theorem 1.}---The maximal maskable set for an arbitrary $d$-dimensional qudit state lays on a hyperdisk in $d^2-1$ dimensional Euclidean space.

	\begin{proof}---
		Let $\rho$ be a density matrix on the Hilbert space $\mathcal{H}_A$, which is the mixed states to be masked, and $\mathcal{U}$ is the masker, a unitary operator acting on $\mathcal{H}_A\otimes\mathcal{H}_B$ as: $\mathcal{U}\rho\otimes|0\rangle\langle0|\mathcal{U}^\dag=\rho^{AB}.$
		Further, let the marginal state $\rho_A=\mathrm{Tr}_B(\rho^{AB})$, then

		\begin{align}
			\rho_A(k,l)=\sum^{d^2-1}_{i=1}a^{(kl)}_ix_i+b^{(kl)}, k,l=1,2,...,d.
		\end{align}

		where $\rho_A(k,l)$ is the $k$-th row and $l$-th column element of matrix $\rho_A$, $a^{(kl)}_i,b^{(kl)}\in \mathbb{C}$.
		A given subset of quantum states in $\mathcal{H}_A$ is maskable iff. $\rho_A(k,l)$ for $\rho_A$, and similarly $\rho_B(k,l)$ for $\rho_B$, are constant, namely, $\exists~ h^{(kl)}_X\in \mathbb{C},~ \rho_X(k,l)\equiv h^{(kl)}_X$, $X=A~ \mathrm{or}~ B$.
		Obviously, for a masker $\mathcal{U}$, if all $\rho_X(k,l)$ are constant for all $\rho$, then masking of all quantum qudit pure states is possible, which contradicts the known results, so $a^{(kl)}_i,i=1,2,...,d$ can't all be zero. In other words, there exists nontrivial hyperplane $\mathrm{Re}\,\rho_X(k,l)=\mathrm{Re}\,h^{(kl)}_X$ or $\mathrm{Im}\,\rho_X(k,l)=\mathrm{Im}\,h^{(kl)}_X$, which contains all of the maskable states corresponding to the masker $\mathcal{U}$. So the maximal maskable set of states with respect to $\mathcal{U}$ always lays on a hyperdisk in $d^2-1$ dimensional space.

	\end{proof}

	Note that for the qubits, we have $d = 2, r_d \equiv1$. In this case, the map between the set of points in the Bloch sphere and the set of qubit states is a bijection, so the maximal maskable set of qubit states with respect to $\mathcal{U}$ is on a disk in the Bloch sphere. Moreover, an arbitrary disk passing through the point $\rho_0=(x_0,y_0,z_0)$ in the Bloch sphere can always be represented as $\mathcal{D}^\theta_\alpha(\rho_0)$.
	We now prove that any qubit disk $\mathcal{D}^\theta_\alpha(\rho_0)$ can be masked. Let's define the qubit masker $\mathcal{U}^\theta_\alpha$, with
	\begin{align}
		\label{eqs:uat}
		\mathcal{U}_\alpha^\theta &= \begin{pmatrix}\cos(\alpha/2) & 0 & e^{-i\theta}\sin(\alpha/2) & 0 \\ 0 & \cos(\alpha/2) & 0 & e^{-i\theta}\sin(\alpha/2) \\ 0 & \sin(\alpha/2) & 0 & -e^{-i\theta}\cos(\alpha/2) \\ \sin(\alpha/2) & 0 & -e^{-i\theta}\cos(\alpha/2) & 0\end{pmatrix}.
	\end{align}
	We apply the masker on the initial state to find the final state resulted by masking, which reads:
	\begin{align}
		\rho^{AB} = \mathcal{U}^\theta_\alpha\rho\otimes|0\rangle\langle0|{\mathcal{U}^\theta_\alpha}^\dag.
	\end{align}
	Then, direct calculation shows that the marginal density matrices for any $\rho\in\mathcal{D}_\alpha^\theta(\rho_0)$ are identical:
	\begin{subequations}
		\begin{align}
			\rho_A\equiv \mathrm{Tr}^B[\rho_{AB}] &= \frac{1}{2}(1+c)|0\rangle\langle0|+\frac{1}{2}(1-c)|1\rangle\langle 1|, \\
			\rho_B\equiv \mathrm{Tr}^A[\rho_{AB}] &=  \frac{1}{2}(1+c)|0\rangle\langle0|+\frac{1}{2}(1-c)|1\rangle\langle 1|, \\
			\text{with}~~~ c &= x_0\sin\alpha\cos\theta+y_0\sin\alpha\sin\theta+z_0\cos\alpha. \nonumber
		\end{align}
	\end{subequations}

	Furthermore, we comment on the measure of the maskable set. Because the volume of the Bloch sphere is $\frac{4\pi}{3}$ and the maximal maskable set is a disk with a volume of zero, the measure of the maskable set corresponding to all qubit states is zero. This conclusion also holds for qudit states, because the volume of all qudit states is greater than the volume of the ball of with $r = \sqrt{\frac{2}{d(d-1)}}$, and the volume of hyperdisk is zero because it is restricted in an $d^2-1$ dimensional space. This lead to the conclusion that a maskable set corresponding to any isometry has zero measure.

	\subsection{Strict inclusion of broadcastable set in maskable set}

	In this section, we consider the relationship between the set of the maskable and broadcastable qudit states. Here, we show that the latter set is a proper subset of the former one. Our start point is that the density matrices of the broadcastable states commute\,\cite{Barnum96}, and the proof is based on the following observation.

	\textit{Lemma 1.}---Commuting mixed states set can be simultaneously masked.

	\begin{proof}
		The density operators of mixed states are represented by Hermitian matrices. As such, the density operators of commuting mixed states can be simultaneously diagonalized. Suppose the spectrum decomposition of $\rho^A_s$ is:

		\begin{align}
			\label{eq:bcable}
			\rho^A_s=\sum^{d}_{i=1}p_i|i\rangle\langle i|, ~\mathrm{for~ all}~ \sum^{d}_{i=1}p_i=1.
		\end{align}

		The equation (\ref{eq:bcable}) represents a hyperplane in $d$-dimensional Euclidean space, and have zero volume measure. Let's construct a Vandermonde matrix $(a_{kl})_{d\times d}$ of unit roots, namely

		\begin{eqnarray}
			(a_{kl})_{d\times d}=\left(\begin{array}{ccccc}
				1 & 1 & 1 & \cdots& 1\\
				x_1 & x_2 & x_3 & \cdots& x_d\\
				x^2_1 & x^2_2 & x^2_3 & \cdots& x^2_d\\
				\vdots & \vdots & \vdots & \ddots& \vdots\\
				x^{d-1}_1 & x^{d-1}_2 & x^{d-1}_3 & \cdots& x^{d-1}_d\\ \end{array} \right
			),
		\end{eqnarray}

		where $x_i$ are the roots of $x^d = 1$. One may verify that the inner product of the two different row vectors of this matrix $(a_{kl})_{d\times d}$ is zero. Let we denote $|\Psi_k\rangle=\frac{1}{\sqrt{d}}\sum^{d}_{l=1} a_{kl}|l\rangle|l\rangle$, then ${\{|\Psi_k\rangle\}^d_{k=1}}$ is a set of orthonormal bases.
		Also, we choose the isometric linear operator $\mathcal{M}$ acting on the base $\{|k\rangle\}^d_{k=1}$ in the following way,
		\begin{eqnarray}
			\label{eqs:vdm}
			& |k\rangle \rightarrow |k\rangle|1\rangle\rightarrow \Psi_k,~ k=1,\ldots,d.\nonumber
		\end{eqnarray}

		From the definition of the isometry, $\mathcal{M}$ can always be dilated to an unitary operator $\mathcal{U}$ on the space $\mathcal{H}_A\otimes \mathcal{H}_B $, that is, there exist a unitary operator $\mathcal{U}$,  $\mathcal{U}|k\rangle|0\rangle=\mathcal{M}|k\rangle=|\Psi_k\rangle$.  Note that $\mathcal{U}\rho^A_s\otimes|0\rangle\langle0|\mathcal{U}^\dag=\rho^{AB}_s$. Direct calculation of $\rho_A$ and $\rho_A$ yields the marginal density matrices
		\begin{eqnarray}\rho_A=\rho_B=\frac{1}{d} \sum^{d}_{k=1}\sum^{d}_{i=1}p_i \parallel a_{ik}\parallel|k\rangle\langle k|\equiv \frac{1}{d}\sum^{d}_{k=1}|k\rangle\langle k|.\nonumber\end{eqnarray}
		Where $\rho_A=\rho_B$ is guaranteed by the symmetric form of the isometry. We see that the eigenvalues $p_i$ cease to appear in the final reduced density matrices, so the masking succeed.
	\end{proof}

	We now proceed to prove the Theorem 2 in the main text.

	\textit{Theorem 2.}---Any qudit broadcastable set is a proper subset of some qudit maskable sets. Moreover, a qudit maskable set has nonzero measure in $d$-dimensional Euclidean space.

	\begin{proof}
		Let
		\begin{align}
			\ket{\psi_0} = \left(\frac{d-1}{d}+\frac{1}{d}i\right)\ket{1} + \sum_{k=2}^d\left(-\frac{1}{d}+\frac{1}{d}i\right)\ket{k},
		\end{align}
		where $i$ is the unit imaginary number. Applying the isometry $\mathcal{M}$ defined in (\ref{eqs:vdm}) on $\ket{\psi_0}$, we find $\mathrm{Tr}_A(\mathcal{M}\ket{\psi_0}\bra{\psi_0}\mathcal{M}^\dagger) = \mathrm{Tr}_B(\mathcal{M}\ket{\psi_0}\bra{\psi_0}\mathcal{M}^\dagger) = \frac{1}{d}\sum^{d}_{k=1}|k\rangle\langle k|$, and thus $\ket{\psi_0}\bra{\psi_0}$ and $\rho_s^A$ is simultaneously maskable by $\mathcal{M}$.
		By the linearity of the trace operation, the set of quantum states living in $d+1$-dimensional real space $\mathbf{P}=\{\rho\,\vert\,\lambda_0 \ket{\psi_0}\bra{\psi_0} + \sum_{i=1}^d \lambda_k \ket{k}\bra{k}, ~\sum_{k=0}^d \lambda_k=1\}$ is maskable by $\mathcal{M}$. It is straightforward to check that for $\prod_{i=0}^d\lambda_i\neq0$ the states $\rho\in\mathbf{P}$ do not mutually commute, so they are predominantly not broadcastable.
		Moreover, as the qudit maskable set lives on the $d+1$-dimensional hyperplane which has nonzero volume measure on the $d$-dimensional Euclidean space, so the dimension of the set of the maskable quantum states is at least 1 plus which of the broadcastable ones.

	\end{proof}

	The Theorem 2 is most intuitive when inspected in the qubit case. The geometric form of a set of commuting mixed states of qubit states in the Bloch sphere is a line segment passing through the center of the sphere. However, any disk in the Bloch sphere is a maskable set, so any broadcastable set is a proper subset of some maskable sets. As the result, the resources of the maskable set are far more abundant than those of the broadcastable set.

	Finally, we comment on the existence of high dimensional maskable set. We have shown in the Theorem 1 that the maximally maskable set lays on a hyperdisk in $d^2-1$-dimensional Euclidean space, that is, the dimension of a maskable set is possible to be as high as $d^2-2$. Generally speaking, the construction of such a set and its corresponding masker is complicated, and remains an open question in the research area of quantum information masking. Here, we takes a step forward by constructing a qudit isometry in even dimension $d=2n$, which masks a $d^2-d$-dimensional set. Consider a group of mappings from $\mathcal{H}^A$ to $\mathcal{H}^A \times \mathcal{H}^B$, acting on the qudit computational bases, $\{\ket{0}, \ket{1}, \ldots, \ket{d-1}\}$, as:
	\begin{align}
		\ket{0} \to \ket{0}^A\ket{0}^B &\to \hspace{29pt} \cos\frac{\alpha_1}{2} \ket{0}^A\ket{0}^B + \sin\frac{\alpha_1}{2} \ket{1}^A\ket{1}^B, \nonumber\\
		\ket{1} \to \ket{1}^A\ket{0}^B &\to e^{-i\theta_1}\left(\sin\frac{\alpha_1}{2} \ket{0}^A\ket{0}^B - \cos\frac{\alpha_1}{2} \ket{1}^A\ket{1}^B\right), \nonumber\\
		\ket{2} \to \ket{2}^A\ket{0}^B &\to \hspace{29pt} \cos\frac{\alpha_2}{2} \ket{2}^A\ket{2}^B + \sin\frac{\alpha_2}{2} \ket{3}^A\ket{3}^B, \nonumber\\
		\ket{3} \to \ket{3}^A\ket{0}^B &\to e^{-i\theta_2}\left(\sin\frac{\alpha_2}{2} \ket{2}^A\ket{2}^B - \cos\frac{\alpha_2}{2} \ket{3}^A\ket{3}^B\right), \nonumber\\
		&\vdots \label{eq:hdmasker}\\
		\ket{d-2} \to \ket{d-2}^A\ket{0}^B &\to \hspace{29pt} \cos\frac{\alpha_n}{2} \ket{d-2}^A\ket{d-2}^B + \sin\frac{\alpha_n}{2} \ket{d-1}^A\ket{d-1}^B, \nonumber\\
		\ket{d-1} \to \ket{d-1}^A\ket{0}^B &\to e^{-i\theta_n}\left(\sin\frac{\alpha_n}{2} \ket{d-2}^A\ket{d-2}^B - \cos\frac{\alpha_n}{2} \ket{d-1}^A\ket{d-1}^B\right). \nonumber
	\end{align}
	The above mappings can be dilated to an isometry. Furthermore, let $\mathcal{D}=\{D_{i}|1\leqslant i\leqslant n\},\, \mathcal{F}=\{F_{jk}|1\leqslant j<k\leqslant n\}$, we can verify that for two groups of fixed real numbers $\{p_1, p_2, \ldots, p_n\}, \{c_1, c_2, \ldots, c_n\}$ with $p_i \geq0, \sum_{i=1}^n p_i=1$, the isometry always masks the following quantum states:
	\begin{align}
		\rho(\mathcal{D}, \mathcal{F}) &= \begin{pmatrix}
			p_1 \mathbf{D}_1 & \mathbf{F}_{12} & \mathbf{F}_{13} & \cdots & \mathbf{F}_{1n} \\
			\mathbf{F}_{12}^\dagger & p_2 \mathbf{D}_2 & \mathbf{F}_{23} & \cdots & \mathbf{F}_{2n} \\
			\mathbf{F}_{13}^\dagger & \mathbf{F}_{23}^\dagger & p_3 \mathbf{D}_3 & \cdots & \mathbf{F}_{3n} \\
			\vdots & \vdots & \vdots & \ddots & \vdots \\
			\mathbf{F}_{1n}^\dagger & \mathbf{F}_{2n}^\dagger & \mathbf{F}_{3n}^\dagger & \cdots & p_n \mathbf{D}_n
		\end{pmatrix},
	\end{align}
	where the $2\times2$ matrices on the diagonal are defined as $\mathbf{D}_i=\frac{1}{2} (I_2+x_i \sigma_x +y_i \sigma_y +z_i \sigma_z)$, with $x_i \sin\alpha_i \cos\theta_i + y_i \sin\alpha_i \sin\theta_i + z_i \cos\alpha_i =c_i$, and the choices for the elements of the $2\times2$ off-diagonal matrices $F_{jk}$ are completely arbitrary, provided that they result in a valid density matrix. To see this more clearly, we directly calculate the marginal states after masking, which are found to be:
	\begin{align}
		\rho_A=\rho_B &=\frac{1}{2}\begin{pmatrix}
			p_1(I_2+c_1\sigma_z) & 0 & 0 & \cdots & 0 \\
			0 & p_2(I_2+c_2\sigma_z) & 0 & \cdots & 0 \\
			0 & 0 & p_3(I_2+c_3\sigma_z) & \cdots & 0 \\
			\vdots & \vdots & \vdots & \ddots & \vdots \\
			0 & 0 & 0 & \cdots & p_n(I_2+c_n\sigma_z)
		\end{pmatrix}.
	\end{align}
	They only depends on $p_i$ and $c_i$, and is irrelevant of the exact forms of $\mathcal{D}$ and $\mathcal{F}$. So, the dimension of the maskable set under this construction is $2|\mathcal{D}|+8|\mathcal{F}| =2n+8(n-1+1)(n-1)/2 =4n^2-2n =d^2-d$. This completes the proof that there exists qudit maskable set that have nonzero measure on $d^2-d$-dimensional Euclidean space.

	\clearpage

	\section{Experimental implementations}

	In the experimental part of the Supplementary Material, we first give a detailed account for the relation between the photon fusion gate and the photonic masking machine. Next, using this relation, we show how the extensibility of our masking machine, that how it can be exploited to mask qudit states, and how a qubit generated from practical single-photon source can be masked using a coherent state produced by a pulsed laser. Finally, we discuss a protocol and the experimental result of a masking-based secret sharing, where a secret image is distributed among three recipients. The reconstruction uses no classical communication between the sender and the receivers, but only classical communication between the receivers about the marginal states of their received qubit.

	\subsection{The relation between the fusion gate and the masking machine}

	Here, we give a detailed analysis of the relation between the fusion gate\,\cite{Browne05} and the masking operator $\mathcal{U}_0^0$ in the picture of second quantization. The schematic illustration of the photon fusion gate is shown in Fig.\,\ref{figs:fuse}, where the wavefunction of two photons interfere on a polarizing beam splitter (PBS).
	The effective surface of the PBS is denoted by the anti-diagonal line crossing the center of the PBS. Regardless of the original direction of propagattion, a photon will be reflected by the PBS if it has vertical polarization (denoted by $\ket{V}$), and will pass through the PBS if its direction of polarization is horizontal (denoted by $\ket{H}$).

	\begin{figure}[htbp]
		\centering
		\includegraphics[width = .25 \textwidth]{fuse.pdf}
		\caption{Schematic illustration of the photon fusion gate.}
		\label{figs:fuse}
	\end{figure}

	A general polarization state $\ket{\psi^A}$ of a single photon can be spanned on the horizontal-vertical basis, namely, $\ket{\psi^A}=\beta \ket{H} + \gamma e^{i\phi} \ket{V}$, with $\beta, \gamma \text{~and~} \phi$ being real numbers and $\beta^2+\gamma^2=1$. For simplicity, in this section we only consider the case of pure states--the situation for mixed states can be proved immediately using the linearity of quantum mechanics.
	In the masking machine (cf. Fig.\,1 in the main text), this photon will interfere with another photon with fixed polarization of $\ket{D}=(\ket{H}+\ket{V})/\sqrt{2}$ on the surface of the PBS. After the PBS gate, the wavefunction propagating from down to is sent to Alice, and the one propagating from left to right is sent to Bob. Moreover, only the photon detection events that result in one photon detected in Alice's site (A) and one in Bob's site (B) are registered as coincidence counting, and the other events are discarded.

	To address the quantum interference on the PBS, we introduce the photon creation operator $a_\mu^{\nu\dagger}$ of the photon, which denotes the creation of a photon with polarization of $\mu\in\{H, V\}$ propagating toward the site $\nu\in\{A, B\}$.
	The conversion of the creation operators of the input photons on the fusion gate reads\,\cite{yhshih}:
	\begin{align}
		\label{eq:modconv}
		a_H^{A\dagger} \to a_H^{A\dagger}, ~~~a_V^{A\dagger} \to ia_V^{B\dagger}, ~~~a_H^{B\dagger} \to a_H^{B\dagger}, ~~~a_V^{B\dagger} \to ia_V^{A\dagger}.
	\end{align}

	Using the creation operators, the two-photon input state can be expressed as:
	\begin{align}
		\ket{\psi^A} \otimes \ket{D} &= (\beta \ket{H} + \gamma e^{i\phi} \ket{V})\otimes \ket{D} = (\beta a_H^{A\dagger} + i\gamma e^{i\phi} a_V^{A\dagger}) \,\frac{1}{\sqrt{2}} (a_H^{B\dagger} + a_V^{B\dagger}) \ket{0^A0^B},
	\end{align}
	where $\ket{0^A 0^B}$ denotes the vacuum state. By substituting (\ref{eq:modconv}) into the above equation, we acquire the form of the two-photon state after the fusion gate:
	\begin{align}
		\ket{\psi^A} \otimes \ket{D} \overset{\text{Fusion}}{\longrightarrow} \ket{\psi^{AB}} &= (\beta a_H^{A\dagger} + i\gamma e^{i\phi} a_V^{B\dagger}) \frac{1}{\sqrt{2}} (a_H^{B\dagger} + i a_V^{A\dagger}) \ket{0^A 0^B} \nonumber \\
		&= \left(\frac{\beta a_H^{A\dagger} a_H^{B\dagger} -\gamma e^{i\phi} a_V^{A\dagger} a_V^{B\dagger}}{\sqrt{2}} + i\frac{\beta a_H^{A\dagger} a_V^{A\dagger} + \gamma e^{i\phi} a_H^{B\dagger} a_V^{B\dagger}}{\sqrt{2}}\right) \ket{0^A 0^B}.
	\end{align}
	The second fraction in the parenthesis, however, does not result in coincidence counting. Consequently, the effective wavefunction after the PBS that reaches coincidence counting, after normalization, reads:
	\begin{align}
		\ket{\psi^{AB}}=\beta \ket{H^A H^B} - \gamma e^{i\phi} \ket{V^A V^B}.
	\end{align}
	Observe that the fusion gate effectively implements an entangling projection $\ket{HH}\bra{HH}-\ket{VV}\bra{VV}$ on the input two-photon state. This conclusion is first drawn in \cite{lmduan06}, and here the form of the projector is slightly different due to the extra $\pi/2$ phase the photons picked up upon reflection\,\cite{yhshih}.

	By checking the two marginal density operators of the bipartite system, we find,
	\begin{align}
		\mathrm{Tr}_A\left(\ket{\psi^{AB}}\bra{\psi^{AB}}\right) &= \mathrm{Tr}_B\left(\ket{\psi^{AB}}\bra{\psi^{AB}}\right) = \begin{pmatrix}
			\beta^2 & 0 \\
			0 & \gamma^2
		\end{pmatrix}.
		\label{eq:dmreduced}
	\end{align}
	We see that the parameter $\phi$ has ceased to appear in (\ref{eq:dmreduced}), which means that the class of quantum states with the same $\beta$, $\gamma$ and different phases $\phi$ are masked. These states lays on a latitudinal parallel on the Bloch sphere. Following the definition of the masker, we conclude that the fusion gate realizes $\mathcal{U}_0^0$.

	Furthermore, we explicitly give the form of the masking isometry $\mathcal{U}$ from the form of mode conversion taking place in the fusion gate. Because the fusion gate and the Hadamard gate acting on the auxiliary photon link the input and the output states as $\begin{pmatrix}\beta \\ \gamma\end{pmatrix} \otimes \begin{pmatrix}1 \\ 0\end{pmatrix} \overset{\text{Fusion}}{\longrightarrow} \begin{pmatrix}\beta \\ 0 \\ 0 \\ -\gamma\end{pmatrix}$, up to a renormalization of the final state, we can choose
	\begin{align}
		\mathcal{U}_0^0=\begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & -1 \\ 0 & 0 & -1 & 0\end{pmatrix}.
	\end{align}
	Note that the from of the second and the fourth columns of the $\mathcal{U}_0^0$ can be freely chosen, as long as the value assignments preserve unitarity of the evolution. Finally, because
	that a $SU(2)$ rotation $\mathcal{R}_\alpha^\theta = \begin{pmatrix}\cos (\alpha/2) & e^{-i\theta}\sin (\alpha/2) \\ \sin (\alpha/2) & -e^{-i\theta}\cos (\alpha/2) \end{pmatrix}$ rotates the edge of the disk $\mathcal{D}_\alpha^\theta (\rho_0)$ to a latitudinal circle on Bloch sphere, we can construct the masking isometry $\mathcal{U}_\alpha^\theta$ as is given in (\ref{eqs:uat}) by applying the $SU(2)$ rotation before the masking gate. By doing this, the final form of the mapping becomes:
	\begin{align}

		\mathcal{U} &= \mathrm{CNOT}\,(\mathcal{R}(\alpha, \theta)\otimes I_2) = \begin{pmatrix}\cos(\alpha/2) & 0 & e^{-i\theta}\sin(\alpha/2) & 0 \\ 0 & \cos(\alpha/2) & 0 & e^{-i\theta}\sin(\alpha/2) \\ 0 & \sin(\alpha/2) & 0 & -e^{-i\theta}\cos(\alpha/2) \\ \sin(\alpha/2) & 0 & -e^{-i\theta}\cos(\alpha/2) & 0\end{pmatrix},
	\end{align}
	which is in accord with the masking isometry defined in (\ref{eqs:uat}).

	\subsection{Extensibility of the masking protocol for qudit}

	The polarization degree of freedom offers a natural photonic qubit state suitable for quantum information processing. In order to realize the masking isometry in higher dimensions, we can introduce more photon pairs to synthesize a composite qudit from some qubits. The quantum circuit for masking higher-dimensional states is shown in Fig.\,\ref{fig:qudit}. Below, we explain the scheme in detail.

	\begin{figure}[htbp]
		\centering
		\includegraphics[width=.4 \textwidth]{qudit.pdf}
		\caption{Applying the masking machine on qudit states. The tilde line denotes the fusion gate.}
		\label{fig:qudit}
	\end{figure}

	A composite qudit $\ket{\psi}=c_0\ket{0}+c_1e^{i\phi_1}\ket{1} +\cdots+c_{d-1}e^{i\phi_{d-1}}\ket{d-1}$ with the dimension of $d\leq2^n$ can be encoded digit-wise on the horizontal ($H$) and vertical ($V$) polarizations of $n$ photons $\ket{\cdot}^1 \ket{\cdot}^2 \cdots \ket{\cdot}^n$, and the correspondence between the qudit and the joint states of the photon polarizations reads
	\begin{align*}
		\ket{0}&\leftrightarrow \ket{HH\ldots HH}, \\
		\ket{1}&\leftrightarrow \ket{HH\ldots HV}, \\
		\ket{2}&\leftrightarrow \ket{HH\ldots VH}, \\
		\ket{3}&\leftrightarrow \ket{HH\ldots VV}, \\
		&\cdots \\
		\ket{2^n-1}&\leftrightarrow \ket{VV\ldots VV}.
	\end{align*}
	Henceforth, for simplicity, we have omitted the superscript of photon indexing. To mask the target qudit we need an auxiliary qudit, again composed of $n$ photons, and set at the polarization state $\ket{D}^{\otimes n}$ with $\ket{D}=(\ket{H}+\ket{V})/\sqrt{2}$. Each of the carrier photon then undergoes a fusion operation with an auxiliary photon on a PBS.
	We write down the initial state of the two-qudit system before the fusion gate on the polarization bases, namely,
	\begin{align}
		\ket{\psi}\otimes\ket{D}^{\otimes n} = \left(\sum_{k=0}^{d-1} c_k e^{i\phi_k} \ket{k}\right)\ket{D}^{\otimes n}.
	\end{align}

	The photon fusion gate effectively conducts a postselection so that one and only one photon emerges from each of the output modes. Consequently, for each fusion gate, the two output photons necessarily have the same polarization. In this case, because the inputs of each fusion gate are still single photon, we simply project the corresponding two-photon states on the $\ket{HH}\bra{HH}-\ket{VV}\bra{VV}$ basis (the minus sign before the $\ket{VV}\bra{VV}$ terms is due to the extra $\pi/2$ phase the photons picked up upon reflection) to acquire the evolutions of the calculation bases, which read:
	\begin{align}
		\ket{0}\otimes\ket{D}^{\otimes n} &=\ket{HH\ldots HH}\left\vert DD\ldots DD\right\rangle \nonumber\\
		&\overset{\text{Fusion}}{\longrightarrow} \phantom{(-1)\times} 2^{-n/2}\ket{HH\ldots HH}\ket{HH\ldots HH} = 2^{-n/2}\ket{0}\ket{0}, \nonumber\\
		\ket{1}\otimes\ket{D}^{\otimes n} &=\ket{HH\ldots HV}\left\vert DD\ldots DD\right\rangle \nonumber\\
		&\overset{\text{Fusion}}{\longrightarrow} (-1)\times 2^{-n/2}\ket{HH\ldots HV}\ket{HH\ldots HV} = 2^{-n/2}\ket{1}\ket{1}, \nonumber\\
		\ket{2}\otimes\ket{D}^{\otimes n} &=\ket{HH\ldots VH}\left\vert DD\ldots DD\right\rangle \nonumber\\
		&\overset{\text{Fusion}}{\longrightarrow} (-1)\times 2^{-n/2}\ket{HH\ldots VH}\ket{HH\ldots VH} = 2^{-n/2}\ket{2}\ket{2}, \nonumber\\
		\ket{3}\otimes\ket{D}^{\otimes n} &=\ket{HH\ldots VV}\left\vert DD\ldots DD\right\rangle \nonumber\\
		&\overset{\text{Fusion}}{\longrightarrow} \phantom{(-1)\times} 2^{-n/2}\ket{HH\ldots VV}\ket{HH\ldots VV} = 2^{-n/2}\ket{3}\ket{3}, \nonumber\\
		&\vdots \nonumber\\
		\ket{k}\otimes\ket{D}^{\otimes n} &=\ket{k}\left\vert DD\ldots DD\right\rangle \overset{\text{Fusion}}{\longrightarrow} (-1)^{\mathcal{P}(k)}\times 2^{-n/2}\ket{k}\ket{k}, \\
		&\vdots \nonumber\\
		\ket{d}\otimes\ket{D}^{\otimes n} &=\ket{d}\left\vert DD\ldots DD\right\rangle \overset{\text{Fusion}}{\longrightarrow} (-1)^{\mathcal{P}(d)}\times 2^{-n/2}\ket{d}\ket{d}. \nonumber
	\end{align}
	Here, we have defined $(-1)^{\mathcal{P}(k)}$ as the parity of $k$. In other words, $\mathcal{P}(k)$ is the count of 1's in the binary form of $k$. The result is based on the observation, that every vertically-polarized photon encoding the computational basis $\ket{k}$ of the composite qudit contributes a $\pi$-phase when affected by a fusion gate. Consequently, the quantum state $\ket{\psi_f}$ resulted by fusing pairwise the carrier and auxiliary photons reads:
	\begin{align}
		\ket{\psi}\otimes\ket{D}^{\otimes n} \overset{\text{Fusion}}{\longrightarrow} \ket{\psi_f} = \sum_{k=0}^{d-1} (-1)^{\mathcal{P}(k)} c_k e^{i\phi_k} \ket{k}\ket{k}.
	\end{align}
	Further, the marginal state for Alice, who receives only one of the output composite qudit, reads:
	\begin{align*}
		\rho_f^\text{A} &= \mathrm{Tr}^{B}[\ket{\psi_f}\bra{\psi_f}] = \sum_{k=0}^{d-1} |c_k|^2 \ket{k}\bra{k},
	\end{align*}
	and, by symmetry, the marginal state for Bob is $\rho_f^\text{B} = \rho_f^\text{A}$.

	We see that the relative phases $e^{i\phi_k}$ in the initial state $\ket{\psi}$ have no observable effect on the two receiving party's marginal state. By definition of quantum information masking, we have masked a proper subset of a qudit Hilbert space. In conclusion, by introducing additional photons, quantum informatiom masking can be extended to arbitrary high dimension using our setup, and the masking isometry acts as $\ket{\psi}=\sum_k c_k \ket{k} \to \ket{\psi_f}=\sum_k (-1)^{\mathcal{P}(k)}c_k \ket{kk}$, with $\mathcal{P}(k)$ being the total number of 1's in the   binary representation of $k$.

	\subsection{Extensibility of the masking protocol for different photon sources}

	In the main article, we demonstrated the photonic quantum information masking using the parametric photons from pumping the nonlinear crystals. However, the scope of the masking machine is not limited by the specific photon source. Since the fusion gate only requires overlapping and indistinguishability of the two-photon wavefunction on the PBS to implement the mode conversion, the masking machine may also accept other types of photon sources as the input. Here, we show that the weak coherent light (e.g., attenuated laser beam) can be exploited to mask the quantum information carried by single photons.

	\begin{figure}[htbp]
		\centering
		\includegraphics[width=.28 \textwidth]{fusec.pdf}
		\caption{Schematic illustration of the fusion gate with a single-photon source and a coherent state as the inputs.}
		\label{fig:fusec}
	\end{figure}

	The schematic illustration of the scenario is shown in Fig.\,\ref{fig:fusec}. Specifically, we are interested to mask a carrier photonic qubit, which is generated by a periodically pumped  semi-deterministic single-photon source, using an auxiliary, pulsed laser which have the same spectral profile as the carrier photon, and the same repetition rate, so a single-photon always meets a coherent pulse at the PBS.
	In this case, because the input modes of the fusion gate are no longer single-photons, it is not possible to directly project the input states on $\ket{HH}\bra{HH}-\ket{VV}\bra{VV}$ to extract the output state, and the situation must be treated in the second quantization picture. We again span the polarization state of the carrier photon on the horizontal-vertical basis, that is, $\ket{\psi^A}=\beta\ket{H}+\gamma e^{i\phi}\ket{V}$, and denote as $p$ the efficiency of the single-photon source to successfully generate a photon. Then, the photon state on the carrier side can be expressed as:
	\begin{align}
		\ket{\psi^A}=\left[\sqrt{1-p}+\sqrt{p}(\beta a_H^{A\dagger}+\gamma e^{i\phi} a_V^{A\dagger})\right]\ket{0^A}.
	\end{align}
	On the other hand, the polarization of photons on the auxiliary side is still set to $\ket{D}$. The second quantization form of the input coherent state reads:
	\begin{align}
		\ket{\alpha}=\sum_{k=0}^{+\infty}\frac{1}{k!}\left(\alpha\, a_D^{B\dagger}\right)^k\ket{0^B}=\sum_{k=0}^{+\infty}\frac{\alpha^k}{k!}\left(\frac{a_H^{B\dagger}+a_V^{B\dagger}}{\sqrt{2}}\right)^k\ket{0^B}.
	\end{align}
	Also, that the coherent laser is weak implies $|\alpha|\ll1$, so we make use of the Maclaurin series expansion up to the second order of $\alpha$ to find:
	\begin{align}
		\ket{\alpha} \approx \left[1+\frac{\alpha}{\sqrt{2}}\left(a_H^{B\dagger}+a_V^{B\dagger}\right)+\frac{\alpha^2}{4}\left(a_H^{B\dagger}a_H^{B\dagger}+2a_H^{B\dagger}a_V^{B\dagger}+a_V^{B\dagger}a_V^{B\dagger}\right)\right]\ket{0^B}.
	\end{align}

	We are now in the position to analyze the efficacy of the masking machine for the scenario of single-photon plus coherent laser input. Applying the mode conversion rule (\ref{eq:modconv}) on the joint input state yields
	\begin{align}
		\ket{\psi^A}\otimes\ket{\alpha} =& \left[\sqrt{1-p}+\sqrt{p}(\beta a_H^{A\dagger}+\gamma e^{i\phi} a_V^{A\dagger})\right] \left[1+\frac{\alpha}{\sqrt{2}}\left(a_H^{B\dagger}+a_V^{B\dagger}\right)+\frac{\alpha^2}{4}\left(a_H^{B\dagger}a_H^{B\dagger}+2a_H^{B\dagger}a_V^{B\dagger}+a_V^{B\dagger}a_V^{B\dagger}\right)\right] \ket{0^A0^B}\nonumber\\ \overset{\text{Fusion}}{\longrightarrow}& \left[\sqrt{1-p}+\sqrt{p}(\beta a_H^{A\dagger}+i\gamma e^{i\phi} a_V^{B\dagger})\right] \left[1+\frac{\alpha}{\sqrt{2}}\left(a_H^{B\dagger}+ia_V^{A\dagger}\right)+\frac{\alpha^2}{4}\left(a_H^{B\dagger}a_H^{B\dagger}+2ia_H^{B\dagger}a_V^{A\dagger}-a_V^{A\dagger}a_V^{A\dagger}\right)\right] \ket{0^A0^B}.
	\end{align}
	However, the system subjects to the postselection from coincidence counting, so the final recorded state only contains the terms which have creation operators for both Alice and Bob. By ignoring the vacuum state, expanding the product and discarding the irrelevant terms that does not give coincidence counting, we find:
	\begin{align}
		\ket{\psi^A}\otimes\ket{\alpha} \overset{\text{Fusion}}{\longrightarrow}& \left[\sqrt{\frac{p}{2}}\,\alpha(\beta a_H^{A\dagger}+i\gamma e^{i\phi}a_V^{B\dagger}) \left(a_H^{B\dagger}+ia_V^{A\dagger}\right) + \frac{\alpha^2\sqrt{1-p}}{4} \left(a_H^{B\dagger}a_H^{B\dagger}+2ia_H^{B\dagger}a_V^{A\dagger}-a_V^{A\dagger}a_V^{A\dagger}\right)\right]\ket{0^A0^B} \nonumber\\
		\overset{\text{coin.}}{\longrightarrow}& \left[\sqrt{\frac{p}{2}}\,\alpha(\beta a_H^{A\dagger} a_H^{B\dagger} -\gamma e^{i\phi}a_V^{A\dagger}a_V^{B\dagger}) + \frac{i\alpha^2\sqrt{1-p}}{2}a_H^{B\dagger}a_V^{A\dagger} \right]\ket{0^A0^B}.
	\end{align}
	So, the final state $\ket{\psi^{AB}}$ detected by coincidence counting can be obtained by applying the creation operators on the vacuum state, and switch back to the notation of first quantization. It reads:
	\begin{align}
		\ket{\psi^{AB}} &= (\beta\ket{HH}-\gamma e^{i\phi}\ket{VV}) + i\sqrt{\frac{1-p}{2p}}a\ket{HV},
	\end{align}
	in which the first term is the desired outcome of the masking machine, and the second term is resulted by the multi-photon terms in the coherent state. The normalization is made for the first term. We see that the noise term vanish at the limit $p\gg a$, i.e., when the mean photon numbers of the coherent state per pulse is much less than the efficiency of the single-photon source. The condition is well attainable with current state-of-the-art deterministic single-photon sources, where some realizations with $p>0.7$ have been reported. Our calculation shows that the photonic quantum information masking machine is also exploitable for the broader class of photon sources, and may find its role in future single-photon source-based quantum information and quantum communication scenarios.

	\subsection{Image reconstruction based on qubit masking}

	In the main text, we have proposed a protocol for a tripartite sharing of a secret qubit state. The sender, Alice, distributes a secret quantum state to three recipients, Bobs, who can only use their marginal state to restrict the masked state onto a disk, and have to cooperate together to decode the secret. Precisely, Alice should mask the secret qubit state $\rho_0$ with three different maskers $\mathcal{U}_{(i)}$ and sends one of the qubits from the resulted $\rho_i^{AB}$ to $\text{Bob}_i$, respectively. The Bobs are informed \textit{a priori} the form of the maskers $\mathcal{U}_{(i)}$. However, from the received marginal state, $\text{Bob}_i$ can only restrict the masked state onto a disk. It is when all Bobs cooperate together and compare their result of the marginal state (using only classical communication) that they can unlock the secret. Upon comparison, the three disks will intersect at a single point in the Bloch sphere, which reveals the concealed information $\rho_0$.
	Our protocol also offers some security against eavesdropping providing that the masker and the marginal state are not simultaneously divulged: fabrication of either of the two components in masking may eventuate in reconstruction of a non-physical state that locates outside the Bloch sphere, which exposes the eavesdropper.

	Furthermore, utilizing the homomorphism between the Bloch representation of a qubit and the hue-saturation-luminosity color space, we can use a qubit state to encode a colored pixel, and share a secret picture pixel-wise between three recipients.
	In our protocol of image sharing, the color of each pixel is determined by three float numbers interpreted by the receivers, ``Bobs''. For this course, a correspondence between the string of numbers decoded from Bobs' quantum states and the reconstructed color has to be established. Here, this is achieved by resorting to the hue-saturation-luminosity (HSL) representation of a color.

	We recall the geometric representation for any mixed qubit state in the Bloch sphere, that is,
	$$\rho(\vec{r}) = \rho(x, y, z) = \frac{1}{2} \left(I_2+x\sigma_x+y\sigma_y+z\sigma_z\right). $$

	For every point in the Bloch sphere, we define the three color parameters, viz., hue, saturation and luminosity, as
	\begin{align*}
		h &= \frac{1}{2} + \frac{\arctan(x, y)}{2\pi}, \\
		s &= \frac{x^2+y^2}{1-z^2}, \\
		l &= \frac{1+z}{2}.
	\end{align*}
	respectively, where $\arctan(x, y)\in[-\pi, \pi)$ is the arc tangent of $y/x$, taking into account which quadrant the point $(x, y)$ is in. The value of hue is undefined at the points $(x, y, z)$ with $x=y=0$, and the saturation is undefined at $z=\pm1$. These singularities will not cause confusion to the final rendition of colors, which is described using the red, green and blue (RGB) values of the colors. The RGB and HSL values of a color is linked by
	\begin{align*}
		\{r, g&, b\} = \{f(0), f(8), f(4)\}, \\
		\text{with}~~~ f(n) &= l-s\min(l, 1-l) \max[-1, \min(k-3, 9-k, 1)], \\
		\text{and}~~~~ k&=\left(n+\frac{6h}{\pi}\right) \mod 12.
	\end{align*}
	\begin{figure}[htbp]
		\centering
		\includegraphics[width = .24 \textwidth]{4supp.png}
		\caption{The one-on-one correspondence between possible colors and quantum state's location in the Bloch sphere.}
		\label{figs:hsl}
	\end{figure}

	The bijection of the RGB colors and the quantum state's location in the Bloch sphere is illustrated in Fig\,\ref{figs:hsl}. The singularities on the $z$-axis resolve because they represent grayscale colors can be solely described by the luminosity.
	In our experiment, for every pixel to be transmitted, Alice encodes the message $(x_0, y_0, z_0)$ into the quantum state $\rho_0$. By applying the set of maskers $\mathcal{U}_{\pi/2}^0, \mathcal{U}_{\pi/2}^{\pi/2}$ and $\mathcal{U}_0^0$ on $\rho_0$, and sending one of the resulted qubits to each Bob, she effectively prepares the marginal states for the three Bobs to be $\begin{pmatrix}(1+x_0)/2 & 0 \\ 0 & (1-x_0)/2\end{pmatrix}$, $\begin{pmatrix}(1+y_0)/2 & 0 \\ 0 & (1-y_0)/2\end{pmatrix}$ and $\begin{pmatrix}(1+z_0)/2 & 0 \\ 0 & (1-z_0)/2\end{pmatrix}$, respectively.
	From the form of the masker and the marginal state, each of the three Bobs can determine a disk in the Bloch sphere that the original state must belong to, and the three disks lays on three mutually orthogonal planes so they necessarily intercept at one point in the Bloch sphere. The intercepting point in turn represents the color information of the transmitted pixel.

	\begin{figure}
		\centering
		\includegraphics[width=.6 \textwidth]{share.pdf}
		\caption{Tripartite secret sharing of an image by masking of mixed states. (a) From their corresponding reduced states and the masker information informed by Alice, every Bob independently interpret the possible disk that contains the masked states. The original state is pinned at the intersection of the three disks, and the location uniquely determines a color using the hue-saturation-luminosity representation. (b) A pixel art of a parrot to be shared. (c) The reconstructed image from the input of three Bobs.}
		\label{fig:share}
	\end{figure}

	In our experiment, Alice adopts $\mathcal{U}_{(1)}=\mathcal{U}_{\pi/2}^0,\mathcal{U}_{(2)}=\mathcal{U}_{\pi/2}^{\pi/2}\text{~and~}\mathcal{U}_{(3)}=\mathcal{U}_0^0$ to mask the quantum state $\rho_0$ corresponding to the color of each pixel. The extraction protocol is displayed in Fig.\,\ref{fig:share}a. From the bipartite state resulted by applying $\mathcal{U}_{(i)}$, one particle is sent to $\text{Bob}_i$ and one is kept by Alice. Because the reduced density matrix for Bob after applying the maskers on $\rho_0$ are $(I_2+x\sigma_z)/2,(I_2+y\sigma_z)/2,~\text{and}~(I_2+z\sigma_z)/2$, respectively, this operation effectively grants each Bobs access to one of the float numbers, $x_0,y_0,\text{~and~}z_0$, which constitutes the original state. By repeatedly applying the above procedure, the entire image is split and transmitted. The Bobs can then assemble their deduced parameters to recover $\rho_0$, which corresponds to a certain color of a pixel. The reconstructed image in Fig.\,\ref{fig:share}c reasonably resembles the original one (cf. Fig.\,\ref{fig:share}b), and the correlation between the two colored images, averaging over red-green-blue channels, is calculated to be $99.35\%$.

	\let\oldaddcontentsline\addcontentslin
	\renewcommand{\addcontentsline}[3]{