text:
\section{Introduction}

Recently, there has been considerable interest in
representing graphs with as few objects as possible.
The idea behind this objective is to keep the
visual complexity of a drawing low for the observer.
The types of objects that have been used are straight-line segments
\cite{desw-dpgfs-CGTA07, dm-dptfs-CCCG14, hkms-dpgfg-JGAA18,
  kms-eaadfs-JGAA18} and circular arcs \cite{s-dgfa-JGAA15,
  hkms-dpgfg-JGAA18}.

Chaplick et al.~\cite{cflrvw-dgflfp-GD16} considered \emph{covering}
straight-line drawings of graphs by
lines or planes
and defined the following new graph parameters.
Let $1 \leq l < d$, and let $G$ be a graph. The \emph{$l$-dimensional affine cover}
number of~$G$ in~$\mathbb R^d$, denoted by $\rho^l_d(G)$, is defined as the minimum
number of $l$-dimensional planes in~$\mathbb R^d$ such that~$G$ has a
crossing-free straight-line drawing that is contained in the union of
these planes.  The \emph{weak} $l$-dimensional affine cover
number of~$G$ in~$\mathbb R^d$, denoted by $\pi^l_d(G)$, is defined similarly
to $\rho^l_d(G)$, but under the weaker restriction that only the vertices
are contained in the union of the planes.
Clearly,  $\pi^l_d(G)\le\rho^l_d(G)$, and
if $l'\le l$ and $d'\le d$ then~$\pi^l_d(G)\le
\pi^{l'}_{d'}(G)$ and $\rho^l_d(G)\le \rho^{l'}_{d'}(G)$.
It turns out that it suffices to study the parameters~$\rho^1_2$, $\rho^1_3$,
$\rho^2_3$, and $\pi^1_2$, $\pi^1_3$, $\pi^2_3$:

\begin{thm}[Collapse of the Affine Hierarchy \cite{cflrvw-dgflfp-GD16}]
  For any integers $1\le l<3\le d$ and for any graph $G$, it holds
  that $\pi^l_d(G)=\pi^l_3(G)$ and $\rho^l_d(G)=\rho^l_3(G)$.
\end{thm}

Disproving a conjecture of Firman et al.~\cite{flsw-ewlc2-GD18},
Eppstein~\cite{Epp19} constructed planar, cubic, 3-connected,
bipartite graphs on $n$ vertices with $\pi^1_2(G) \geq n^{1/3}$.
Answering a question of Chaplick et al.~\cite{cflrvw-dgflfp-GD16}
he also
constructed a family of subcubic series-parallel graphs with unbounded
$\pi^1_2$-value.  Felsner~\cite{Fe19} proved that, for every
4-connected plane triangulation~$G$ on $n$ vertices, it holds that
$\pi^1_2(G) \leq \sqrt{2n}$.
Chaplick et al.~\cite{cflrvw-cdgfl-WADS17} also investigated the
complexity of computing the affine cover numbers.  Among others, they
showed that in 3D, for $l \in \{1,2\}$, it is NP-complete to decide
whether $\pi^l_3 (G) \leq 2$ for a given graph~$G$.  In 2D, the
question has still been open, but a related question was raised by
Dujmovi\'c et al.~\cite{dpw-tlg-DMTCS04} already in 2004.  They
investigated so-called \emph{track layouts} which are defined as
follows.  A graph admits a $k$-track layout if its vertices can be
partitioned into $k$ ordered independent subsets such that any pair of
subsets induces a plane graph (w.r.t.\ the order of the subsets).  The
track number of a graph~$G$, $\tn(G)$, is the smallest~$k$ such that
$G$ admits a $k$-track layout.  See also \cite{DJMMUW} for some recent
developments.
Note that in general
$\pi^1_2(G) \ne \tn(G)$; for example, $\pi^1_2(K_4)=2$, whereas
$\tn(K_4)=4$.  Note further that a 3-track layout is necessarily plane
(which is not the case for $k$-track layouts with $k>3$).  Dujmovi\'c
posed the computational complexity of $k$-track layout as an open
question.

While it is easy to decide efficiently whether a graph admits a
2-track layout, Bannister et al.~\cite{bddew-tllpd-Algorithmica19}
answered the open question of Dujmovi\'c et al.\ already for 3-track
layouts in the affirmative.  They first showed that a graph has a
leveled planar drawing if and only if it is bipartite and has a
3-track layout.  Combining this results with the NP-hardness of level
planarity, proven by Heath and Rosenberg~\cite{hr-loguq-SICOMP92},
immediately showed that it is NP-hard to decide whether a given graph
has a 3-track layout.  For $k>3$, deciding the existence of a
$k$-track layout is NP-hard, too, since it suffices to add to the
given graph $k-3$ new vertices each of which is incident to all
original vertices of the graph~\cite{bddew-tllpd-Algorithmica19}.

\paragraph{Our contribution.}
We investigate several problems concerning the {\em weak line cover number} $\pi_2^1(G)$
and the {\em strong plane cover number} $\rho_3^2(G)$:
\begin{itemize}
\item We settle the open question of Chaplick et
  al.~\cite[p.~268]{cflrvw-cdgfl-WADS17}
  by showing that
  it is NP-hard to test whether, for a given planar graph~$G$,
  $\pi^1_2(G)=2$; see Section~\ref{sec:NP}.
\item We show

  that $G_d$, the universal stacked triangulation of depth~$d$, (which
  has treewidth~3) has $\pi^1_2(G_d)=d+1=\log_3(2n_d-5)+1$, where
  $n_d$ is the number of vertices of~$G_d$; see
  Section~\ref{sec:threetrees}.
\item Eppstein has identified classes of treewidth-2 graphs with
  unbounded $\pi^1_2$-value.  We give an easy direct argument showing
  that some 2-tree $H_d$ with $n'_d$ vertices has $\pi^1_2(H_d)\in
  \Omega(\log n'_d)$; see \lncsarxiv{the full
    version~\cite{arxiv}}{Appendix~\ref{sec:2trees}}.
\item Concerning 3D, we show that any $n$-vertex graph~$G$ with
  $\rho^2_3(G)=2$ has at most $5n-19$ edges; see
  Section~\ref{sec:twoplanes}.  This bound is tight.
\end{itemize}

\section{Complexity of Computing Weak Line Covers in 2D}
\label{sec:NP}

In this section we investigate the computational complexity of
deciding whether a graph can be drawn on two lines.

\begin{thm}
  It is NP-hard to decide whether a given plane (or planar) graph~$G$
  admits a drawing with $\pi^1_2(G)=2$.
\end{thm}

\begin{proof}
  Our proof is by reduction from the problem \textsc{Level Planarity},
  which Heath and Rosenberg~\cite{hr-loguq-SICOMP92} proved to be
  NP-hard.  The problem is defined as follows.  A planar graph~$G$ is
  \emph{leveled-planar} if its vertex set can be partitioned into sets
  $V_1,\dots,V_m$ such that $G$ has a planar straight-line drawing
  where, for every $i \in \{1,\dots,m\}$, vertices in~$V_i$ lie on the
  vertical line $\ell_i \colon y=i$ and each edge~$v_j v_k$ of~$G$
  connects two vertices on consecutive lines (that is, $|j-k|=1$).

  Chaplick et al.~\cite{cflrvw-dgflfp-GD16} have shown that every
  leveled-planar graph can be drawn on
  two lines.  The converse, however, is not true.  For example, $K_4$
  is not leveled-planar, but $\pi^1_2(K_4)=2$.

  Therefore, we modify the given graph in three ways.  (a)~We replace
  each edge of~$G$ by a $K_{2,4}$-gadget where the two nodes in one
  set of the bipartition replace the endpoints of the former edge; see
  Fig.~\ref{fig:K24-gadget}.  (b)~We add to the graph~$G'$ that
  resulted from the previous step a new subgraph~$G_0$ (two copies
  of~$K_4$ sharing exactly two vertices), which we
  connect by a path to a vertex on the outer face of~$G$.  (If the
  outer face is not fixed, we can try each vertex.)  In
  Fig.~\ref{fig:spiral}, $G_0$ is yellow
  and the path is red.  The
  length~$L$ of the path is any upper bound on the number of levels
  of~$G'$, e.g., the diameter of~$G'$ (plus~1).  (c)~We attach
  to~$G_0$ a triangulated spiral~$S$ (dark green in
  Fig.~\ref{fig:spiral}).  The spiral makes $L+2$ right turns; its
  final vertex is identified with the outermost vertex of the previous
  turn.  Hence, apart from its many triangular faces, the graph
  $S+G_0$ has a large inner face~$F$ of degree $2(L+2)$ and a quadrangular
  outer face.  Let~$G''$ be the resulting graph.
  It remains to show that $G$ is leveled-planar if and only if $\pi_2^1(G'')=2$.

  ``$\Rightarrow$'': Fix a
  leveled-planar drawing of~$G$.

By doubling the layers and using the new layers to place the large sides of $K_{2,4}$'s, one easily sees that $G'$ is also leveled-planar, see Fig.~\ref{fig:K24-gadget}.
  As shown in Fig.~\ref{fig:spiral}, the large inner face~$F$
  of~$S+G_0$ can be drawn so that it partitions the half{}lines emanating from the origin
  into~$L$ levels.  (It is no problem that consecutive levels are
  turned by $90^\circ$.)  Since we chose~$L$ large enough (in
  particular $L \ge 2m-1$), we can easily draw~$G'$ inside~$F$.  Note
  that the red path attached to~$G_0$ is long enough to reach any
  vertex on the outer face of~$G'$.  Hence, $\pi^1_2(G'')=2$.

  \begin{figure}[tb]
    \begin{subfigure}[b]{.44\textwidth}
      \centering
      \includegraphics{K24-gadget}
      \caption{transforming~$G$ to~$G'$}
      \label{fig:K24-gadget}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.5\textwidth}
      \centering
      \includegraphics{spiral}
      \caption{transforming~$G'$ to~$G''$}
      \label{fig:spiral}
    \end{subfigure}
    \caption{Our reduction from \textsc{Level Planarity}}
    \label{fig:reduction}
  \end{figure}

  ``$\Leftarrow$'': Fix a drawing of~$G''$ on two lines.  The two
  lines cannot be parallel since $G''$ contains $K_{2,4}$ and is not
  outer-planar; so after translation and/or skew we may assume that
  these two lines are the two coordinate axes.
  It is not hard to verify that~$G_0$ must be drawn such that the
  origin is in its interior, at the common edge of the two $K_4$'s.
  Furthermore, given this drawing of~$G_0$, the 3-connected spiral~$S$
  must be drawn as in Fig.~\ref{fig:spiral}.  Due to
  planarity and the fact that~$G'$ is connected to~$G_0$ via the red
  path, $G'$ can only be drawn in the interior of~$F$.  The drawing of
  $S+G_0$ partitions the half{}lines emanating from the origin into
  levels, which we number $1,2,\dots$ starting from the innermost
  level that contains a vertex of~$G'$.  Inside this face, the only
  way to draw the $K_{2,4}$-gadgets is as in
  Fig.~\ref{fig:K24-gadget}, spanning three consecutive levels.  This
  forces all vertices of~$G$ to be placed on the odd-numbered levels
  and the vertices in $G'-G$ on the even-numbered levels.  Now we can
  get a level assignment for~$G$ by reverting the transformation in
  Fig.~\ref{fig:K24-gadget}.

  Hence, $G$ is leveled-planar.

  This shows that our reduction is correct.  It runs in polynomial
  time.
\end{proof}

\section{Weak Line Covers of Planar 3-Trees in 2D}
\label{sec:threetrees}

In this section we consider the weak line cover number~$\pi^1_2$ for
planar graphs,
i.e., we are interested in crossing-free straight-line drawings
with vertices located on a small collection of lines.
Clearly $\pi^1_2(G)=1$ if and only if $G$ is a forest of paths.
The set of graphs with $\pi^1_2(G)=2$, however, is already
surprisingly rich, it contains all trees, outerplanar graphs and
subgraphs of grids~\cite{bddew-tllpd-Algorithmica19,Epp19}.

\emph{Stacked triangulations}, a.k.a.~planar 3-trees or Apollonian
networks, are obtained from a triangle by repeatedly selecting a
triangular face~$T$ and adding a new vertex (the {\em vertex stacked
  inside $T$}) inside~$T$ with edges to the vertices of~$T$.  This
subdivides~$T$ into three smaller triangles, the \emph{children}
of~$T$.

For $d\geq 0$ let~$G_d$ be the \emph{universal stacked
  triangulation of depth $d$}, defined as follows.
The graph~$G_0$ is a triangle~$T_0$, and~$G_d$ (for $d\geq 1$) is
obtained from~$G_{d-1}$ by adding a
stack vertex in each bounded face of~$G_{d-1}$.
Graph~$G_d$ has $n_d=\frac{1}{2}(3^d+5)$
vertices and $3^d$ bounded faces.  We show that its weak line cover
number is $d+1=\log_3(2n_d-5)+1 \in \Theta(\log n_d)$.
(A lower bound of $d$
can also be found in Eppstein's recent book~\cite[Thm. 16.13]{EppBook}.)

\begin{thm}\label{thm:stacked}
  For $d\geq 1$ it holds that $\pi(G_d) = d+1$.
\end{thm}

\begin{proof}
Here we prove only the lower bound; the construction for the upper
bound is illustrated in Fig.~\ref{fig:stacked} and given in
\lncsarxiv{the full version}{Appendix~\ref{app:stacked}}.
Let $\cal{L}$ be a family of lines covering the vertices of a drawing
of~$G_d$.  Let~$a$, $b$, and~$c$ be the vertices of~$T_0$.  We first argue
that at least $d$ lines are needed to cover $V\setminus T_0$.
Let~$x_1$ be stacked into~$T_0$.  There is a line
$L_1\in \cal{L}$ covering $x_1$.  Note that $L_1$ can intersect only
two of the three child triangles of~$T_0$
(where ``intersect'' here means ``in the interior''). Let~$T_1$ be a
child triangle avoided by~$L_1$, and let $x_2$ be the vertex stacked
into~$T_1$.  There is a line $L_2\in \cal{L}$ covering $x_2$.
Let~$T_2$ be a child triangle of~$T_1$ avoided by~$L_2$. Iterating this
yields $d$ pairwise distinct lines in~$\cal{L}$.

To find one additional line in~$\cal{L}$, we distinguish some cases.
If a line $L\in \cal{L}$ covers two vertices of~$T_0$, then it covers
no inner vertex, and we are done.

Assume some line $L_a\in \cal{L}$ intersects~$x_1$ and one vertex
of~$T_0$, say $a$.  Let~$L_b$ and~$L_c$ be the lines intersecting~$b$
and~$c$.  The lines~$L_a$, $L_b$, and~$L_c$ are pairwise different,
else we are in the previous case.  Of the three child triangles of
$T_0$, at most one is intersected by $L_a$ and at most two each are
intersected by $L_b$ and~$L_c$.  Therefore, some child triangle~$T_1$
of~$T_0$ is intersected by at most one of~$L_a$, $L_b$, or~$L_c$.  The
graph~$G_{d-1}$ inside~$T_1$ requires at least $d-1$ lines for its
interior points, and at most one of those lines is~$L_a$, $L_b$,
or~$L_c$, so in total at least $d+1$ lines are needed.

The argument is similar if no line covers two of $a$, $b$, $c$,
and~$x_1$.  The four distinct lines supporting $a$, $b$, $c$,
and~$x_1$ then intersect at most two child triangles each.  So one
child triangle~$T_1$ is intersected by at most two of these lines.
Combining the $d-1$ lines needed for the interior of~$T_1$ with the
two lines that do not intersect it, shows that $d+1$ lines are
needed.
\end{proof}

\begin{figure}[tb]
  \centering
  \includegraphics[scale=.18,trim=800 0 0 0,clip]{stacked}
  \caption{A drawing of $G_2$ that can be extended to a drawing of
    $G_3$ on 5 parallel lines.}
  \label{fig:stacked}
\end{figure}

\section{Maximal Graphs on Two Planes in 3D}
\label{sec:twoplanes}

We now switch to dimension $d=3$ and the strong cover number.
Obviously any graph $G$ with a drawing that is covered by two planes
has at most $6n-12$ edges since it is the union of two planar graphs.
Using maximality arguments and
counting, we show that in fact $G$
has at most $5n-19$ edges if $n\geq 7$.
(The restriction $n \ge 7$ is required since for
$n=3,4,5,6$ we can have $3,6,9,12$ edges.)

We argue first that our bound is tight.
The {\em spine} is the intersection of two planes $A$ and $B$.
Put a path with $n-4$ vertices on the spine.  Add
one vertex in each of the four halfplanes and connect each of these
vertices to all vertices on the spine and to the vertex on the
opposite halfplane\lncsarxiv{. (We provide a figure in the full
  version.)}{; see Fig.~\ref{fig:tightness} in
  Appendix~\ref{app:boundarycases}.}
This yields $n-5$ edges on the path and $2(n-4)+1$ edges in each of
the two planes, so $5n-19$ edges in total.

\begin{thm}
  \label{thm:twoplanes}
  Any graph~$G$ with $\rho^2_3(G)=2$ and $n\geq 7$ vertices has at
  most $5n-19$ edges.
\end{thm}
\begin{proof}
Fix a drawing of~$G$ on planes $A$ and $B$, inducing planar
graphs~$G_A$ and~$G_B$ within those planes.  Let $G_A^+$ and $G_B^+$
be the graphs obtained
from~$G_A$ and~$G_B$ by adding any edge that can be inserted without
crossing, within the same plane, and with at most one bend on the spine.
Clearly it suffices to argue that $G_A^+$ and $G_B^+$ together have at
most $5n-19$ edges.
Let $s$ be the number of vertices on the spine, let $a$ be the number
of vertices of~$G_A^+$ not on the spine, and let $b$ be the
number of vertices of~$G_B^+$ not on the spine.  Clearly, $a+b+s=n$.
We may assume $a\leq b$.
We also assume that $1\leq s\leq n-4$ and that at least one edge
of~$G_A^+$ crosses the spine (so $2\leq a\leq b$); see
\lncsarxiv{the full version}{Appendix~\ref{app:boundarycases}}.

Let $t$ be the number of edges drawn along the spine.  These are the
only edges that belong to $G_A^+$ and $G_B^+$.  Since $G_A^+$ and $G_B^+$
have at least three vertices each, we can bound the number of edges
of~$G$, $m(G)$, as follows:
\begin{align}
\label{equ:1}
m(G) &\le m(G_A^+) + m(G_B^+)-t \,\le\, 3(s+a)-6 + 3(s+b)-6 - t \\
&= 3n-12 + 3s - t \,\le\, 4n-16+2s-t. \nonumber
\end{align}
So we must show that $2s-t\leq n-3$.  Let an {\em internal gap} be a
line segment connecting two consecutive, non-adjacent vertices on the spine.
There are $s-t-1$ internal gaps.  Let the {\em external gap} be the two
infinite parts of the spine.  Note that at least one edge of $G_A^+$ must
cross the external gap, because $G_A^+$ has at least one vertex on each side
of the gap, and we could connect the extreme such vertices (or re-route an existing
edge) to cross the external gap, perhaps using a bend on the spine.
We may further assume that even after such re-routing every internal gap is crossed by at
least one edge of~$G_A^+$.  Otherwise we could delete all edges of~$G_B^+$
passing through the gap, insert the edge between the spine vertices, and
re-triangulate the drawing of $G_B^+$ where we removed edges.  This
would remove an internal gap, but would not decrease the number of edges.
Since no edge can cross two gaps, at least $s-t$ edges of $G_A^+$ cross gaps.
These edges form a planar bipartite graph with
at most $a$ vertices; therefore $s-t\leq 2a-3$.%
\footnote{One might be tempted to
write a bound of $2a-4$ here, but we must allow for the possibility of
$a=2$, in case of which the planar bipartite graph may have $1=2a-3$ edges.}
This yields
$2s-t \leq s + 2a-3 \leq s + a + b - 3 = n-3$
as desired.
\end{proof}

\noindent
We conjecture that the following more general statement
holds:

\begin{quote}
  Any $n$-vertex graph $G$ with $\rho^2_3(G)=k$ has at most
  $(2k+1)(n-2k)+k-1$ edges, for all large enough $n$.
\end{quote}

\paragraph{Acknowledgments.}
This research started at the Bertinoro Workshop on Graph Drawing
2017. We thank the organizers and other participants,
in particular Will Evans, Sylvain Lazard, Pavel Valtr, Sue Whitesides,
and Steve Wismath.  We also thank Alex Pilz and Piotr Micek for
enlightening conversations.

\bibliographystyle{splncs04}


text:
\section{Introduction}

Signatures for the presence of dust are found nearly everywhere in
astrophysics. In the context of star and planet formation, dust is
abundant in molecular clouds and in the circumstellar environments of
a large fraction of stellar objects in the early stages of their
evolution. In the circumstellar disks encircling these young stars,
where planets are thought to form, the interplay between dust and gas
is of paramount importance.

At short wavelengths, dust grains efficiently absorb, scatter, and
polarise the starlight. How much radiation is scattered and absorbed
is a function of both the geometry of the disk and the properties of
the dust. In turns, the amount of absorbed radiation sets the temperature of the
dust (and gas) and defines the amount of radiation that is re-emitted
at longer, thermal, wavelengths.

The last decade has witnessed the improvement of imaging capabilities
with the advent of potent instruments in the optical and
near-infrared, and large millimeter interferometers, providing
detailed views of the disks around young stars. The sensitivity and
wavelength range covered by new instruments is increasing steadily and
the mid- and far-infrared ranges are now being explored efficiently by
the {\sl Spitzer Space Telescope}, and new facilities like Herschel
and ALMA will soon complete the coverage.

With this unprecedented wealth of data, from optical to radio, fine
studies of the dust content and evolution of disks become
possible and powerful radiative transfer (RT) codes are needed to fully
exploit the data.  In this paper we describe such a code, MCFOST.

 That code was used
extensively to produce synthetic images of the scattered light from
disks around young stars. Examples include the circumbinary ring of
GG~Tau \citep{McCabe02, Duchene04}, the large silhouette disk
associated with IRAS 04158+2805 \citep{Menard05}, and an analysis of
the circular polarisation in GSS~30 \citep{Chrysostomo97}.

In \S2 below, we briefly describe MCFOST. In section \S3, tests and
validation of MCFOST are presented. Two examples of applications are
presented in \S4. Firstly, a study is presented of the minimum
mass of disks detectable by {\sl Spitzer} around young solar-like stars, young
low-mass stars, and young brown dwarfs. Secondly, the colours of
parametric disks are compared to the {\sl Spitzer}
colours, for both IRAC and MIPS, of samples of T Tauri stars located in
nearby star forming regions.

\section{Description of the numerical code}
MCFOST is a 3D continuum radiative transfer code based on the
Monte-Carlo method. It was originally developed by \cite{Menard} to
model the scattered light in dense dusty media (including linear and
circular polarisations). In this paper, we present an extended version
of the original code that includes dust heating and continuum thermal
re-emission.

\subsection{Geometry of the computation grid}

MCFOST uses a spatial grid that is defined in cylindrical coordinates,
well adapted to the geometry of circumstellar disks. $N_r$
logarithmically-spaced radial grid cells and $N_z$ linearly-spaced
vertical grid cells are used. The size of the vertical cells follows
the flaring of the disk, with a cutoff at 10~times the scale height
for each radius. MCFOST allows for the density distribution to be
defined arbitrarily, in 3 dimensions, with the limitation that within
each cell, quantities are held constant.
The density distribution of the media under consideration can be
described either by using a parametric model (as in
\S~\ref{sec:bench}) or by including any ad hoc density table, calculated
by other means. In the following, the term ``local'' refers to the
properties of a single cell.

\subsection{Position-dependent dust distributions}
An explicit spatial dependence of the size (and/or composition)
distribution $f(a,\vec{r})$ is implemented in MCFOST.  The dust
properties are therefore defined locally, i.e., each cell of the disk
may contain its own independent dust population. This allows for
example to model dust settling towards the disk midplane, variations
of the chemical composition of the dust, and/or an increase of ice
mantles from the inner, hot regions to the outer, cold edge of the
disk.

\subsection{Optical dust properties}
The dust optical properties are currently computed with the Mie
theory, i.e., grains are spherical and homogeneous.  The optical
properties at a given position in the disk are derived in accordance
with the local size and composition distributions $f(a,\vec{r})$. The
extinction and scattering opacities are given by \begin{equation}
\kappa^\mathrm{ext/sca}(\lambda,\vec{r}) =
\int_{a_\mathrm{min}}^{a_\mathrm{max}} \pi a^2
Q_\mathrm{ext/sca}(\lambda,a) f(a,\vec{r}) \,\mathrm{d}a \end{equation} where $
Q_\mathrm{sca}(\lambda,a)$ and $ Q_\mathrm{ext}(\lambda,a)$ are
respectively the scattering and extinction cross sections of a grain
of size $a$ at a wavelength $\lambda$.

\subsection{The radiative transfer scheme in MCFOST}

The Monte Carlo method allows to follow individual monochromatic
``photon packets'' that propagate through the circumstellar
environment until they exit the computation grid.  The propagation
process is governed by scattering, absorption and re-emission events
that are controlled by the optical properties of the medium (opacity,
albedo, scattering phase function, etc...) and by the temperature
distribution. Upon leaving the circumstellar environment, ``photon
packets'' are used to build an SED and/or synthetic images, one image
per wavelength.

\subsubsection{The emission of ``photon packets''}

In MCFOST, two general radiation sources are considered~: photospheric
emission and circumstellar dust thermal emission.

Any number of stars can be considered, at any position in the
calculation grid. The photospheric emission can be represented by (i)
a sphere radiating uniformly or (ii) a limb-darkened sphere or (iii) a
point-like source, whichever is relevant for the problem under
consideration. The spectrum of photospheric photons can follow a
blackbody or any observed or calculated spectrum. Hot and/or cold
spots can be added to the photospheres.

For the dust thermal emission, the density, temperature and opacities
are considered constant within each grid cell. Thermal photons
emission locations are uniformly distributed within each cell and the
emission is assumed to be isotropic.

\subsubsection{Distance between interactions}

It is natural within a Monte Carlo scheme to estimate the distance a
photon packet ``travels'' between two interactions by means of optical
depth (directly related to the density) rather than by physical,
linear distance. From a site of interaction, i.e., scattering and
emission, the optical depth $\tau$ to the next site is randomly chosen
following the probability distribution $p(\tau) = \exp(-\tau)$.

The distance $l$ is computed by integrating the infinitesimal
optical depth $\kappa^\mathrm{ext}(\lambda,\vec{r}) \rho(\vec{r})
\,\mathrm{d}s$ until the following equality is verified~:
\begin{equation}
\tau_\lambda = \int_0^l \kappa^\mathrm{ext}(\lambda,\vec{r}) \rho(\vec{r}) \,\mathrm{d}s~.
\end{equation}

Once the position of interaction $\vec{r}$ is determined, the
probability that the interaction is a scattering event rather than an
absorption event is estimated with the local albedo~:
\begin{equation}
{p_\mathrm{sca} = \mathcal{A}(\lambda,\vec{r}) = \frac{ \int_{a_\mathrm{min}}^{a_\mathrm{max}} \pi a^2
Q_\mathrm{sca}(\lambda,a) f(a,\vec{r}) \,\mathrm{d}a}{\int_{a_\mathrm{min}}^{a_\mathrm{max}} \pi a^2
Q_\mathrm{ext}(\lambda,a) f(a,\vec{r}) \,\mathrm{d}a}}~.
\end{equation}

\subsubsection{Scattering}

MCFOST includes a complete treatment of linear and circular
polarisations by using the Stokes formalism. The state of a light
packet is described by its Stokes vector, with $I$ representing the intensity, $Q$
and $U$ the linear polarisation, and $V$ the circular
polarisation.  The interaction of a photon with a dust particle is
described by a $4\times4$ matrix, the Mueller matrix, ${\bf S}$, also
known as the scattering matrix. Since the matrix is defined
  with respect to the scattering plane, appropriate rotation matrices
  are required to transform the Stokes parameters between multiple
  scattering events. In the simplifying case of Mie
scattering
, the matrix becomes block-diagonal with only 4 different
non-zero elements, see Eq.~\ref{eq:mueller}.

\begin{equation}
\label{eq:mueller}
\left( \begin{array}{r}
I \\
Q \\
U \\
V \\
\end{array} \right)_{i+1}
=
\left( \begin{array}{rrrr}
S_{11} & S_{12} & 0 & 0 \\
S_{12} & S_{11} & 0 & 0 \\
0    &  0  & S_{33} & S_{34} \\
0 & 0 & -S_{34} & S_{33} \\
\end{array} \right)
\left( \begin{array}{r}
I \\
Q \\
U \\
V \\
\end{array} \right)_{i}
\end{equation}

The direction of scattering is defined by two angles, the scattering
angle $\theta$ and the azimuth angle $\phi$ in spherical
coordinates. Each individual element $S_{ij}$ of the matrix is a
function of the scattering angle and of the wavelength, $S_{ij} =
S_{ij}(\theta,\lambda)$.  The scattering angle $\theta$ is randomly
chosen from the pretabulated cumulative distribution function~: \begin{equation}
\label{eq:phase_theta}
F(\theta) = \frac{\int_0^\theta S_{11} (\theta')\,\sin \theta' \mathrm{d}\theta'}{\int_0^\pi S_{11} (\theta')\,\sin \theta' \mathrm{d}\theta'}.
\end{equation}

For linearly unpolarised incoming light packets, i.e., $Q,U = 0$, the
distribution of azimuthal angle $\phi$ is isotropic. For light packets
with a non-zero linear polarisation, $P =\sqrt{Q^2+U^2}/I$, the
azimuthal angle is defined relative to its direction of polarisation
and determined by means of the following cumulative distribution function~:
\begin{equation}
\label{eq:phase_phi}
F_\theta(\phi) =
\frac{1}{2\pi}\left(\phi-\frac{S_{11}(\theta)-S_{12}(\theta)}{S_{11}(\theta)+S_{12}(\theta)}\,P\,\frac{\sin(2\phi)}{2}\right)
\end{equation}
where $\theta$ was previously chosen from Eq.~\ref{eq:phase_theta}
\citep{Solc89}.

The Mueller matrix used in Eqs. \ref{eq:mueller},
\ref{eq:phase_theta}, and  \ref{eq:phase_phi} can be defined in two different
ways : one Mueller matrix per grain size or one mean Mueller matrix per
grid cell. In the latter case, the local Mueller matrix
$S(\lambda,\vec{r})$ is defined by~:
\begin{equation}
\label{eq:def_S}
{\bf S}(\lambda,\vec{r}) = \int_{a_\mathrm{min}}^{a_\mathrm{max}}
{\bf S}(\lambda,a) f(a,\vec{r}) \,\mathrm{d}a
\end{equation}
where $S(\lambda,a)$ is the Mueller matrix of a grain of size $a$ at a wavelength $\lambda$.

In the former case, the grain size must be
chosen explicitly for each event following the probability law~:
\begin{equation}
p(a)\mathrm{d} a = \frac{\pi a^2
  Q_{\mathrm{sca}}(a)f(a,\vec{r})\,\mathrm{d} a}{\int_{a_\mathrm{min}}^{a_\mathrm{max}}\pi a^2 Q_{\mathrm{sca}}(a)f(a,\vec{r})\,\mathrm{d} a}~.
\end{equation}
The two methods are strictly equivalent. They can be used alternatively
to minimize either the memory space or the computation time required.

\subsubsection{Absorption and radiative equilibrium\label{sec:eq_rad}}
The temperature of the dust particles is determined by assuming
radiative equilibrium in the whole model volume, and by assuming that
the dust opacities do not depend on temperature.

The dust thermal balance should take into account the
thermal coupling between gas and dust. In high density regions, e.g.,
close to the disk midplane, the coupling is very good and the gas
temperature is expected to be close to the dust temperature. At the
surface layers of the disk, on the other hand, densities become much
lower and gas-dust thermal exchanges are reduced.

Two extreme assumptions are implemented in MCFOST~: either the
gas-dust thermal exchange is perfect and gas and dust are in local
thermodynamic equilibrium (LTE) throughout the disk, or there is no
thermal coupling between gas and dust at all. In the latter case, each
grain size has its own, different temperature, despite being in
equilibrium with the same radiation field.

Under the assumption of LTE, all dust grains have the same
temperature, which is also equal to the gas temperature. In this case, the
(unique) temperature of the cell is given by the radiative
equilibrium equation~:
\begin{equation}
4\pi\int_0^\infty \kappa_i^\mathrm{abs}(\lambda) B_\lambda(T_i)
\,\mathrm{d}\lambda = \Gamma_i^\mathrm{abs}
\end{equation}
where $\Gamma_i^\mathrm{abs}$ is the  energy absorption rate. If only
passive heating is considered, this rate
can be written from the mean intensity, leading to~:
\begin{equation}
\label{eq:radequil}
4\pi\int_0^\infty \kappa_i^\mathrm{abs}(\lambda) B_\lambda(T_i)
\,\mathrm{d}\lambda = 4\pi\int_0^\infty \kappa_i^\mathrm{abs}(\lambda)
J_\lambda \,\mathrm{d}\lambda.
\end{equation}
Any extra internal source of energy (viscous heating for instance) can
be easily included by adding the corresponding term on the
right hand side of Eq.~\ref{eq:radequil}.

Each time a packet $\gamma$ of a given wavelength $\lambda$ travels
through a cell, the quantity $\Delta l_\gamma$, the distance
travelled by the packet in the cell, is computed. The mean intensity
$J_\lambda$ in the cell is then derived following \citet{Lucy99}:
\begin{equation}
J_\lambda = \frac{1}{4\pi V_i} \sum_{\gamma} \varepsilon\, \Delta
l_\gamma = \frac{1}{4\pi V_i} \sum_{\gamma} \frac{L_*}{N_\gamma}\,
\Delta l_\gamma
\end{equation}
where $V_i$ is the volume of cell $i$ and $\varepsilon$ the individual
luminosity of a packet. The radiative equilibrium equation can then be
written as~:
\begin{equation}
\int_0^\infty \kappa_i^\mathrm{abs}(\lambda)
B_\lambda(T_i) \, \mathrm{d}\lambda = \frac{L_*}{4\pi V_i N_\gamma} \sum_{\lambda,
\gamma} \kappa_i^\mathrm{abs}(\lambda) \Delta l_\gamma~.
\label{eq:Lucy}
\end{equation}
This method is known as the {\sl mean intensity} method.

When gas and
dust, as well as dust grains of different sizes, are not thermally
coupled,
 there is a different temperature for each grain size and the
radiative equilibrium must be written independently for
each of them~:
\begin{equation}
4\pi\int_0^\infty \kappa_{i}^\mathrm{abs}(\lambda,a) B_\lambda(T_{i}(a))
\,\mathrm{d}\lambda = \Gamma_i^\mathrm{abs}(a)
\end{equation}
where $\kappa_{i}^\mathrm{abs}(\lambda,a)$, $T_{i}(a)$ and $
\Gamma_i^\mathrm{abs}(a)$ are the
opacity, temperature and energy absorption rate of the dust grains of size $a$ in cell $i$.

The calculation of $\int_0^\infty \kappa_i^\mathrm{abs}(\lambda)
B_\lambda(T_i) \, \mathrm{d}\lambda = \sigma T^4 \kappa_P / \pi$ ($\kappa_P$
is the Planck mean opacity) is very time consuming and we pretabulate
these values at $N_T = 1\,000$ logarithmically spaced temperatures
ranging from $1$K to the sublimation temperature, which is of order
$1\,500$ K. The temperature $T_i$ of each cell is obtained by interpolation
between the pretabulated temperatures.

To speed up calculations, \citet{Bjorkman01} introduced the concepts
of {\sl immediate reemission} and {\sl temperature correction}. They
are implemented in MCFOST.  In essence, when a packet is absorbed, it
is immediately re-emitted and its wavelength chosen by taking into
account the temperature correction given by the following probability
distribution~:
\begin{equation} p_\lambda \, \mathrm{d}\lambda \propto
\kappa_\mathrm{i}^\mathrm{abs}(\lambda)
\left(\frac{\mathrm{d}B_\lambda}{\mathrm{d}T}\right)_{T_i} \, \mathrm{d}\lambda~.
\label{eq:correct_PDF}
\end{equation}

\citet{Baes05} suggested that this wavelength distribution adjustment
method may not be entirely safe when used in combination with the
{\sl mean intensity} writing of the absorption rate described above.
Instead, MCFOST determines the cell temperatures $T_i$ used in
Eq. \ref{eq:correct_PDF} by solving~:
\begin{equation}
\int_0^\infty \kappa_i^\mathrm{abs}(\lambda)
B_\lambda(T_i) \, \mathrm{d}\lambda = \frac{L_*}{4\pi V_i N_\gamma} N_{\gamma\
  \mathrm{abs}, i}
\label{eq:Lucy_alt}
\end{equation}
where $N_{\gamma\ \mathrm{abs}}$ is the current number of
absorbed photons in cell $i$. Eq. \ref{eq:Lucy} is used only at the
end,
 when all photons have been propagated, to improve the accuracy of
the temperature determination.  Using Eq.~\ref{eq:Lucy_alt} instead of
Eq.~\ref{eq:Lucy} when computing temperature for
Eq.~\ref{eq:correct_PDF} is in no way a limitation because, at this
stage of the computation, the aim is to compute the radiation field
and not the whole temperature structure. Of course, the temperature in
optically thin cells will not be precisely determined because little
absorption occurs in these cells. On the other hand, as a consequence,
little reemission will also occur and the contribution of these
reemission events will remain very small. By construction of the
\cite{Bjorkman01} method, temperature is best converged in the cells
that contribute most to the radiation field. So, when all photons have
been propagated, temperature corrections are no longer needed and the
temperature of all cells can be computed following the \cite{Lucy99}
method, providing a much higher final accuracy.

\subsubsection{SED calculation: sampling of the radiation field\label{sec:sampling}}

To improve its computational efficiency when producing SEDs, MCFOST
uses different samplings of the radiation field, i.e., different
strategies to set the energy of a photon packet.  The energy of a
packet, hence the number of photons it contains, is chosen and
optimised depending on the goals of the calculation. On the one hand,
the convergence of the temperature distribution is optimized when all
photon packets have the same {\sl energy}, independently of their
wavelengths. This procedure insures that more photon packets are
emitted in the more luminous bins. On the other hand, the computation
of SEDs itself is more efficient when it is the {\sl number} of photon
packets that is held constant for all wavelength bins instead.  In
that case it is the energy of the packets that is wavelength
dependent.  Thus, MCFOST is made more efficient by computing the
temperature and SEDs with a two-step process:
\begin{itemize}
\item {\bf Step 1} is the temperature determination. Photon packets
are generated and calculated one by one at the
stellar photosphere and propagated until they exit the computation
grid, i.e., until they exit the circumstellar medium.  Upon scattering,
the propagation vector of a packet is modified, but not its
wavelength. Upon absorption however, packets are immediately
re-emitted, in situ and isotropically, but at a different wavelength
calculated according to the temperature of the grid cell. For
this re-emission process, the concept of {\sl immediate reemission}
and the associated {\sl temperature correction method} proposed by
\cite{Bjorkman01} are used. In this step, all photon packets have the
same luminosity
  \begin{equation}
    \label{eq:E_step1}
    \epsilon = L_*/N_{\gamma step1}
  \end{equation}
where $L_*$ is the bolometric luminosity of the star and $ N_{\gamma
step1}$ the number of packets generated,
 and are randomly scattered/absorbed within the disk. The
concept of {\sl mean intensity} suggested by \cite{Lucy99} is further
used to compute the final temperature at the end of the step. This
allows to efficiently reduce noise in the temperature estimation for
optically thin cases. Step 1 allows for a fast convergence of the
temperature but is rather time-consuming when used to derive SEDs,
especially in the low-energy, long wavelength regime. To speed up
convergence, a different scheme (step) is used to compute SEDs.

\item {\bf Step 2} computes the SED  from the
final temperature distribution calculated in step 1.  In step 2, the number
of photon packets $N_{\gamma
step2}$ is held constant at all wavelengths.
Step 2 therefore maintains a comparable noise level in each wavelength
bins and significantly reduces convergence time by limiting the CPU
time spent in high luminosity bins and focusing on low luminosity
bins. In this step, photon packets are always scattered, no absorption
is allowed. However, at each interaction, the Stokes vector is
weighted by the probability of scattering, $p_\mathrm{sca}$, to take
into account the energy that would have been removed by
absorption. This allows all the photons to exit the disk and to
contribute to the SED, although with a reduced weight.
Packets at a given wavelength $\lambda$ are randomly
emitted from the star and the disk, by cells $i$, with the respective
probabilities~:
\begin{equation}
p_* = \frac{L_*(\lambda)}{L_*(\lambda)~+~\sum_i~w_i(\lambda)\,L_i(\lambda)}~,
\end{equation}
\begin{equation}
{\rm and}~ p_i = \frac{w_i(\lambda)\,L_i(\lambda)}{L_*(\lambda)~+~\sum_i~w_i(\lambda)\,L_i(\lambda)}~.
\end{equation}
where the luminosities are defined by~:
\begin{equation}
L_*(\lambda) = 4 \pi^2R_*^2\,B_\lambda(T_*)~,
\end{equation}
\begin{equation}
{\rm and}~ L_i(\lambda) = 4 \pi \, m_i\, \kappa_i^{\mathrm{abs}}(\lambda) \,B_\lambda(T_i)~.
\end{equation}
In order to avoid generating photon packets so deep in the (very
optically thick) disk that
none of them would escape with an appreciable energy left,  a ``dark
zone'' is defined for each wavelength. This dark zone includes all
cells from which the optical depth is at least $\tau_\lambda = 30$, in any
direction, to get out of the computationnal grid\footnote{The value
  $\tau_\lambda = 30$ was found to be a good compromise between
  result accuracy and computational time.}. No
photon packet is emitted from this zone, the emission probability from cells
inside the dark zone being set to zero by fixing a weight
$w_i(\lambda)=0$. A weight $w_i(\lambda) = 1$ is given to cells above the dark zone
limit. Furthermore, packets that would enter
the dark zone during their random walk are killed, because they have
no chance to exit the model volume with significant energy.
The luminosity of the $N_{\gamma
step2}$ packets emitted at a given wavelength $\lambda$ is determined
by the total energy that the star and the disk radiates at this
wavelength~:
\begin{equation}
\epsilon_\lambda = \frac{L_*(\lambda)~+~\sum_i~w_i(\lambda)\,L_i(\lambda)}
{N_{\gamma step2}}~.
\end{equation}
\end{itemize}

\subsubsection{Intensity and polarisation maps}

A scheme similar to step 2 used above to calculate SEDs is deployed
in MCFOST to calculate synthetic maps. Maps are calculated, one
wavelength at a time. Scattering is treated explicitly for all
wavelengths but again, absorption is replaced by the appropriate
weighting of the intensity to avoid losses of photon packets.  Images
are produced, in all four Stokes parameters, by classifying photon
packets that exit the calculation grid into inclination bins,
themselves divided into image pixels. The same process applies to
scattered light images in the optical and thermal emission maps in the
far infrared and radio range. The relative contribution of the photospheric and dust thermal
emission is dependent on the wavelength.  It is for instance critical
 in the
  mid-infrared and beyond where it is the (scattered)
thermal emission from the inner disk that starts to dominate over the
photospheric contribution. The correct fraction is calculated
explicitly through the temperature and density of each disk cell and
each contribution can be mapped separately for
comparison. Interferometric visibilities
and phases can be obtained by Fourier transform of images.

\begin{figure}[t]
\centering
\includegraphics[width=8cm,angle=270]{T_r.eps}
\caption{Plots of the radial temperature (upper panel) and the
difference (lower panel) for the most optically thick case,
$\tau_\mathrm{V} = 100$, arbitrarily using RADICAL as the
reference.
In both panels, the results of MCFOST are represented by the thick
solid line. Thin solid lines are the results from MC3D, dot-dashed
lines from MCTRANSF, dotted lines from RADMC and dashed lines from
STEINRAY. In the upper panel, because all curves are very similar,
MCFOST has been shifted by 200K for clarity.
\label{fig:Temp_radiale}}
\end{figure}

\section{\label{sec:test} Tests and validation of MCFOST}

Numerous tests have been performed to check the efficiency, stability,
and accuracy of MCFOST. In this section, we present a comparison of
MCFOST with the benchmark results published by \cite{Pascucci04}
(hereafter P04).

\begin{figure}[t]
\includegraphics[angle=270,width=8.5cm]{diff_sed100.eps}
\caption{\label{fig:diff} Plots of the difference between the model
SEDs for the most optically thick case, $\tau_\mathrm{V} = 100$,
arbitrarily using RADICAL as the reference code. The linestyle-to-code
pairing is the same as in Fig.~\ref{fig:Temp_radiale}.
Three different inclination angles are considered: $i=12.5^\circ$
(upper panel), $i=42.5^\circ$ (middle panel), and $i=77.5^\circ$
(lower panel).
}
\end{figure}

\subsection{\label{sec:bench} Benchmark definition}
 The geometry tested involves a central point-like source radiating as
a T=$5\,800$K blackbody with a bolometric luminosity $L = L_\odot$
encircled by a disk of well defined geometry and dust content. The
disk extends from 1 AU to 1\,000 AU. It includes spherical dust grains
made of astronomical silicate. Grains have a single size radius of
$0.12\ \mu$m and a density of $3.6\ \mathrm{g.cm}^{-3}$. The optical
data are taken from \cite{Draine84}. The disk is flared with a
Gaussian vertical profile. Radial power-law distributions are used for
the surface density and the scale height (see P04 for more details).

To compare with P04, the scattering is also forced to be isotropic and
polarisation is not calculated, i.e, all the scattering information is
contained in the scattering cross section alone, $Q_{sca}$, acting on
the I Stokes parameter only.

The number of grid cells is set to $N_r = 50$ and $N_z = 20$ in the
radial and vertical directions, respectively. We used $10^7$ photon
packets to calculate the temperature distribution (step 1) and
$2\,10^6$ photon packets per wavelength to generate the SED (step
2). The total run-time is $20$ minutes for the most optically thick
case on a bi-processor (Intel Xeon) computer running at a clock rate of
2.4 GHz. The run-time memory space needed is $10$ MBytes.

\subsection{Results}

The results calculated with MCFOST for all the test cases presented in
P04 are available at
{\url{http://www-laog.obs.ujf-grenoble.fr/~cpinte/mcfost/}}. Only the
most optically thick case ($\tau_\mathrm{V}=100$, see P04), i.e., the
most difficult one, is presented in Figs.~\ref{fig:Temp_radiale} and
\ref{fig:diff}.

The temperature profiles in the disk and shape of the emerging SEDs
calculated with MCFOST are in excellent agreement with the calculation
presented in P04. Differences in the radial temperatures computed by
MCFOST do not exceed 4\% with respect to the other 5 codes in the
optically thin case. For the most optically thick case, differences
are generally of the order of 2-5\% (see Fig.~\ref{fig:Temp_radiale}).
The results from MCFOST are always within the range of those produced
by other codes.

Regarding SEDs, MCFOST produces results comparable to all other codes
to better than 10\% for low optical thickness. Deviations as large as
15\% are however observed in the most defavorable case, i.e., high
tilt and optical thickness (either at very short wavelength or in the
10 $\mu$m silicate feature),
but again the results lie in the range of the values produced by
other codes (see Fig.~\ref{fig:diff}).

Various tests were further performed to show the independence of the
MCFOST results with respect to grid sampling, inclination
sampling, and position of the vertical density cut-off in the disk.
Results are not noticeably affected. A complete description of all
numerical tests is presented in \cite{Greta05}.

\subsection{Energy conservation at high optical depth}
As discussed in section \ref{sec:sampling}, SEDs can be constructed from packets escaping the system in both steps
of the computation (the SED from step 2 being
better converged by construction). As photon packets are generated
independently through different processes in step 1 and 2, our adaptive sampling of the
radiation field method allows an \emph{a posteriori} check of
the accuracy of calculations.

By integrating over all inclination
angles and wavelengths,
the emerging flux must correspond in both cases to the star luminosity
(in case of passive heating). This is true by construction for the SED
calculated in step 1.

Assuming a standard flared disk geometry and ISM like grains, and
varying (via the disk mass) the V band equatorial optical depth from $0.1$ to $10^{9}$, we
find that the luminosities calculated through the 2 methods always
agree to within better than $1\%$.

\section{Example of applications of MCFOST}

In this section we present calculations of the SEDs and infrared
photometric colours of young objects surrounded by circumstellar
disks. The colours calculated by MCFOST are compared with photometric
observations obtained with the {\sl Spitzer Space
Telescope}. Calculated infrared fluxes as a function of disk mass are
also compared with {\sl Spitzer}'s detection limits, as obtained for
example during legacy surveys like {\sl c2d} \citep{Evans03}.  The goals are (i) to explore
whether a disk with ``representative'' geometrical parameters can
match the observed infrared colours of solar-like young stars and (ii)
also to estimate the range of disk masses detectable by {\sl Spitzer}
in such programmes.

\subsection{Disk geometry and central object properties}

\begin{table}[t]
  \centering
  \caption{Star parameters}
  \begin{tabular}{l|ccc}
    Central object & Mass ($M_\odot$) & Radius ($R_\odot$) &
    $T_\mathrm{eff}$\\
    \hline
    T Tauri  star & 1.0 & 1.94 & 4\,420 K\\
    Low mass star  & 0.2 & 1.28 & 3\,365 K\\
    Brown dwarf & 0.05  & 0.57 & 2\,846 K\\
    \hline
  \end{tabular}
  \label{tab:parameters}
\end{table}

\begin{figure*}[t]
  \includegraphics[angle=270,width=0.325\textwidth]{sed_tt.eps}
  \hfill
  \includegraphics[angle=270,width=0.325\textwidth]{sed_vlms.eps}
  \hfill
  \includegraphics[angle=270,width=0.325\textwidth]{sed_bd.eps}
\caption{\label{fig:SED_Spitzer} SEDs for different dust disk masses
  and central objects, for a pole-on disk. The left panel presents the
  results for the T Tauri star and $R_\mathrm{in}  = 0.1$  AU, the central panel for the very low
  mass star and the right panel for the brown dwarf (defined in table \ref{tab:parameters}). The line type
  corresponds to the disk dust mass (in solar masses). The horizontal
  bars indicate {\sl Spitzer}'s detection limits in the four IRAC
  bands and  the $24\ \mu$m  and  $71\ \mu$m MIPS bands
  for the molecular cloud mapping programme of the {\sl Cores to
  Disks} legacy survey. Additional MIPS limits (light gray bars~:
  $170$ mJy at $24\ \mu$m and $1\,000$ mJy at $71\ \mu$m) for deeper
  pointed observations of the programme are also plotted at $24\ \mu$m
  and $71\ \mu$m.}
\end{figure*}

We assume a disk geometry that is representative of the geometries
derived from images of several disk sources (HK Tau B~:
\citealp{Stapelfeld98} ; HV Tau C~: \citealp{Stapelfeldt03} ; IRAS 04158+2805~:
\citealp{Glauser06} ; TW Hya~: \citealp{Krist00} ; LkHa 263 C~:
\citealp{Chauvin02}). The disk geometry is kept constant through the
various model runs but the disk mass is varied. Only disks with \emph{dust
masses} smaller than $10^{-6} M_\odot$ are considered in order to focus
on the transition between optically thick and optically thin disks,
{\sl i.e.}, from class II to class III young stellar objects.

The disk is defined by a flared geometry with a vertical Gaussian
profile $\rho(r,z) = \rho_0(r)\,\exp(-z^2/2\,h(r)^2)$.  Power-law
distributions describe the surface density,
$\Sigma(r)=\Sigma_0\,(r/r_0)^{\alpha}$, and the scale height, $ h(r) =
h_0\,(r/r_0)^{\beta}$, where $r$ is the radial coordinate in the
equatorial plane, $h_0 = 10$ AU is the scale height at the reference
radius $r_0 = 100$ AU, $\beta=1.125$, and $\alpha=-1$. The
outer radius is set to $R_\mathrm{out} = 300$ AU.  Two values of
the disk inner radius are considered: $0.1$~AU and $1$~AU.
A grid of models is calculated from dust disk masses ranging from
$10^{-13}$ to $10^{-6} M_\odot$. Calculations are presented with
masses sampled logarithmically, by steps of $\sqrt{10}$.

A grain size distribution $\mathrm{d}n(a) \propto a^{-3.7}
\,\mathrm{d}a$ of spherical particles ranging from $a_\mathrm{min} =
0.03\ \mu$m to $a_\mathrm{max} = 1\ $mm, and optical constants from
\cite{Mathis89}, namely their model A, are used.

 The central stars are assumed to emit like blackbodies.
Three
different central light sources are considered. They are chosen to
represent a typical T Tauri star, a low mass star, and a brown dwarf
(see Tab.~\ref{tab:parameters} for parameters). The radii and
effective temperatures of the T Tauri and low-mass stars are obtained
from the models of \cite{Siess00}, assuming an age of $2$ Myrs and a
metallicity $Z = 0.1$. For the brown dwarf, a mass of $0.05$
M$_\odot$ is selected. The photospheric properties
are taken from \cite{Chabrier00} and from \cite{Baraffe02} dusty
models. These parameters are summarized in Tab.~\ref{tab:parameters}.

All
models are calculated for a distance d=140 pc, typical of nearby star
forming regions.

\subsection{Detectability of disk infrared excesses with {\sl Spitzer}}

In Fig.~\ref{fig:SED_Spitzer}, SEDs are presented as a function of
disk dust mass for the three objects to evaluate the minimum
detectable disk mass as a function of the luminosity of the central
star. For comparison, a set of detection limits reached by {\sl Spitzer} for
large mapping programmes is overplotted. Here we take the values
quoted by the {\sl Spitzer} legacy survey team {\sl Cores to Disks} for
their molecular cloud mapping programme \citep{Evans03} \footnote{See the {\sl c2d} website for
more information \url{http://peggysue.as.utexas.edu/SIRTF/}}.

At all 4 IRAC bands, the photospheric emission is always above the
detection limit, for all types of objects. At $24\ \mu$m, the T
Tauri and low-mass photospheres are detected. However, the photosphere
of the brown dwarf is $\approx 3$ times lower than the detection
limit. This is of course dependent on the exact detection limit chosen
and on the actual mass of the brown dwarf.

For the case considered here, i.e., detection limits for large mapping
programmes and substellar object of 0.05M$_\odot$, a flux detection at
$24\ \mu$m readily implies the presence of circumstellar material. A
more careful analysis is required to confirm a disk excess in the case
of more massive objects. At $71\ \mu$m, the photosphere is detected
for none of the objects and the presence of circumstellar dust can be
directly inferred from any detection.

\begin{figure*}[th]
  \includegraphics[width=8cm]{c1.eps} \hfill
  \includegraphics[width=8cm]{c2.eps} \caption{IRAC colour-colour diagrams. Black filled
  circles corresponds to models with $R_\mathrm{in} = 0.1$ AU and
  large grains ($a_\mathrm{max} = 1$ mm), gray triangles to models
  with $R_\mathrm{in} = 1.0$ AU and large grains, white squares to
  models with $R_\mathrm{in} = 0.1$ AU, with interstellar dust
  ($a_\mathrm{max} = 1\ \mu$m). Large symbols show inclinations
  smaller than $80^\circ$ and small symbols $i > 80^\circ$. Class II
  and Class III data points from \cite{Hartmann05} were converted to
  surface densities using a Gaussian convolution ($\sigma = 0.07$
  mag). The solid and dashed contours define the location of class II
  and class III objects; respectively. The levels are [0.1, 0.2, 0.4,
  0.8] of the maximum of density.}  \label{fig:IRAC_colours}
\end{figure*}

Reverting the argument, for an observing programme with such detection
limits, {\sl Spitzer} is expected to detect, at $71~\mu$m, dust disk masses as low as
$3\,10^{-9}$ M$_\odot$ around T Tauri stars, with $R_\mathrm{in}=0.1$ AU (Fig. \ref{fig:SED_Spitzer}, left panel). For low mass stars, disk
masses greater than $3\,10^{-8}$ M$_\odot$ will be detected (Fig. \ref{fig:SED_Spitzer}, central panel), and for a
brown dwarf of 0.05M$\odot$, disk masses larger than $3\,10^{-7}$
M$_\odot$ are required for
detection (Fig. \ref{fig:SED_Spitzer}, right panel). All results are given in terms of {\sl dust mass}. A
gas-to-dust ratio has to be assumed to extrapolate to total disk
masses.

\subsection{Spitzer colours}
\cite{Hartmann05} have used {\sl Spitzer}/IRAC to determine the
location of TTS (Class II and III) in colour-colour diagrams. Here we
attempt to reproduce the global location of these sources in the
colour-colour diagrams. The precise location depends on the model
geometry.  Similar studies of IRAC colours have been presented
elsewhere \citep{D'Alessio06}. Here, we present the colours of a few
parametric models of passive disks to explore the range of parameters
that reproduce IRAC colours. We also compare the same model with MIPS
colours recently published \citep{Padgett06}.

\begin{figure}[th]
  \includegraphics[width=8cm]{c3.eps}
 \caption{MIPS
  colour-colour diagram. The symbols are the same as in
  Fig.~\ref{fig:IRAC_colours}. The solid and dashed contours define the
  location of class II and class III from \cite{Padgett06},
  respectively. The Gaussian convolution was performed with $\sigma =
  0.07$ mag along the H - K axis and $\sigma = 0.35$ mag along the
  K -[24] axis.
  \label{fig:MIPS_colours}}
\end{figure}

Fig. \ref{fig:IRAC_colours} and \ref{fig:MIPS_colours} show the modelled
IRAC and MIPS colours for the previously defined T Tauri ($1 M_\odot$)
models.  Models for central objects with lower masses (not displayed
for clarity) presents similar colours and, as a consequence, the
comparison of our results for a unique central object mass with the
distribution of masses of \cite{Hartmann05} and \cite{Padgett06}
sources is relevant.
Additional models with interstellar dust, where
$a_\mathrm{max} = 1\ \mu$m and larger inner radius $R_\mathrm{in}=1$
AU are also presented. Three data sets ( $R_\mathrm{in}=0.1$ AU and
$a_\mathrm{max} = 1$\,mm ; $R_\mathrm{in}=1$ AU and $a_\mathrm{max} =
1$\,mm ; $R_\mathrm{in}=0.1$ AU and $a_\mathrm{max} = 1\ \mu$m) are
then displayed, each one calculated for $15$ different disk masses and
for $21$ different inclinations, linearly sampled from pole-on to
edge-on.  Contour plots showing the observed colours of Taurus class II
and class III sources (from \citealp{Hartmann05} in
Fig.~\ref{fig:IRAC_colours} and \citealp{Padgett06} in Fig.~\ref{fig:MIPS_colours}) are surimposed.

Broadly speaking, models sample the two regions of the colour-colour
diagram where sources are found~: low disk mass models corresponding
to class III objects and higher disk mass models to class II.

The slight offsets seen between the lowest mass models and the class
III region is due to the colours of the blackbody spectrum we used.
Using a more realistic spectrum would produce, for instance, a [3.6] -
[4.5] colour index of $-0.04$ mag instead of $0.1$ mag for the blackbody spectrum, in better agreement with the observations.

For increasing disk masses, all IRAC and MIPS colour indices increase
toward the class II region. The models also populate the intermediate
region where no object are observed. These results can be used to
confirm that the transition between class II (disk dust mass $>
3\,10^{-9}$ M$_\odot$) and class III (disk dust mass $< 10^{-10}$)
M$_\odot$ must be fast enough so that no or very few objects are
detected in samples of nearly 50 sources \citep{Hartmann05}.

A clear difference is seen between models depending on the disk inner
radius~: models with $R_\mathrm{in}=1$ AU produce redder [4.5] -
[5.8] and [5.8] - [8.0] colours.  In particular, with $R_\mathrm{in}=1$
AU the [5.8] - [8.0] colour index becomes larger than the observed
value by about $0.5$ mag for large masses, suggesting that
$R_\mathrm{in}=1$ AU is an upper limit for most class II disks in
Taurus. This suggestion is supported by the H - K / K - [24]
diagram, where no excess is seen in H - K for models with $R_\mathrm{in}=1$
AU, well below the average H - K $=0.3$ seen in observations
(see gray triangles in Fig.~\ref{fig:IRAC_colours} and \ref{fig:MIPS_colours}).
  Such a behaviour is explained by the much higher
temperature reached in models with $R_\mathrm{in}=0.1$ AU, for which
disk emission peaks between $2$ and $3\, \mu$m, whereas it peaks
around $10\, \mu$m in models with $R_\mathrm{in}=1.0$ AU.

Models close to edge-on ($i > 80^\circ$) tend to produce redder colour
at large masses than their less tilted counterparts. Because the
central star becomes more attenuated by its disk, the contribution
from the disk emission progressively dominates, leading to rising
mid-infrared SEDs. The resulting colours are reminiscent of Class I
objects, and highlight the necessity to determine the geometry and
inclination, in order to correctly assess the nature of a class I
source, {\sl i.e.}, the presence or not of a massive envelope
\citep{Kenyon95,White04}

No significant influence of the grains size is detectable for models
where the star is seen directly (with $i < 80^\circ$).  At high
inclinations, on the contrary, models with interstellar dust produce
redder colour indices than models with large grains, for which
absorption tends to be gray.

For completely edge-on models, when no direct light from the star
reaches the observer, short wavelength colours~: [3.6] - [4.5] and H
- K strongly depend on scattering properties.  Models with
interstellar dust present almost no excess at short wavelength, most
of the light being scattered starlight. At longer wavelength, direct
emission from the disk is seen, producing red indices. For models with
larger grains, the albedo remains non negligible in the infrared, and
the contribution of thermal emission from the central (and hidden)
parts of the disk scattered by the outer parts towards the observer is
non negligible and adds to the stellar scattered light, resulting in
redder [3.6] - [4.5] and H - K indices than with interstellar
dust.

\section{Summary}

We have presented a new continuum 3D radiative transfer code, MCFOST. The efficiency and reliability of MCFOST was tested considering the benchmark configuration defined by P04.
 MCFOST was shown to calculate temperature
distributions and SEDs that are in excellent agreement with previous
results from other codes.

Sets of models for young solar-like stars and low-stars are presented
and compared to the {\sl Spitzer}'s detection limits for the
programme of the {\sl Cores to Disks} legacy survey.  Minimal disks
masses of $\approx 10^{-9}$ M$_\odot$ for a T Tauri star and $\approx
10^{-7}$ M$_\odot$ for a brown dwarf are needed for the disk to be
detectable by {\sl Spitzer} in the mapping mode.

IRAC and MIPS colours of Taurus Class II and III objects are compared
to passively heated models with a ``representative'' disk
geometry. The average location of the two classes of objects are well
reproduced, as well as the extreme colours of some of the objects, that
may correspond to highly tilted disks.  $R_\mathrm{in} =1$ AU is found
to be maximum inner radius required to account for the observed colours
of Class II objects. Further modelling of individual sources,
combining optical, infrared and millimeter photometry with images
and/or IRS spectroscopy, will be needed to better understand
individual disks, their
dust properties and their evolution.

\begin{acknowledgements}

Computations presented in this paper were performed at the Service
Commun de Calcul Intensif de l'Observatoire de Grenoble (SCCI). We thank the {\sl Programme National de Physique Stellaire}
(PNPS) and {\sl l'Action Sp\'ecifique en Simulations Num\'eriques pour
l'Astronomie} (ASSNA) of CNRS/INSU, France, for supporting part of
this research. Finally, we wish to thank the referee, C.P. Dullemond, for
his comments which have helped to improve the manuscript.

\end{acknowledgements}

\bibliographystyle{aa}

redpajama_set_name: RedPajamaArXiv


text:
\section*{Introduction}

Efficient coding of quantum information relies on the capability to generate high-quality states of light with non-classical photon statistics, both in the continuous-~\cite{Lvovsky2020} and discrete-variable~\cite{Zhong2020} regimes. In particular, states with negative Wigner functions have shown to be crucial for implementing sophisticated protocols and information processing, such as error-correction~\cite{Gottesman2001} and one-way quantum computing~\cite{Zheng2021}. Different strategies for generating such states were reported, relying for instance on high-order nonlinearity~\cite{Rempe2019,Kosuke2022,FLY22}, qubit interaction~\cite{Hastrup2021} and linear optical interaction ~\cite{Lund2004,Etesse2014,Weigand2018,Eaton2019}. A class of protocols of particular interest is that of iterative schemes \cite{Fiurasek2005,Etesse14,Etesse2014,Lund2004,Weigand2018,Eaton2019}, where photonic `resource' states are combined together in an iterative fashion, allowing for growing step-by-step the size of the states. However, the experimentally implemented protocols are to date restricted to few operations, due to protocol scalability issues arising from the simultaneous need for multiple heralding signals and from the increase of experimental overhead~\cite{Etesse2015,Sychev17}. This ultimately limits the average photon number in the generated state.\\
In this article, we propose a novel scheme to overcome the current experimental limitations, by relying on a single shot operation with a spectrally multiplexed resource. Our proposal takes advantage of the multimode aspect of state-of-the-art experimentally investigated quantum photonic sources in order to greatly enhance the capabilities of iterative generation protocols. This improvement occurs at two levels: (i) the number of operations to implement is reduced to one, and, as a consequence, (ii) the success probability of the overall operation is strongly increased. To illustrate our proposal, we show that Schr\"odinger cat state (SCS) and GKP-codes \cite{Gottesman2001} can be generated in an efficient manner, proving the validity of the approach, and opening a new path in the domain of non-Gaussian photonic state preparation.\\
In Sec.~\ref{sec_prot} we recall the principle of iterative protocols and propose a fully parallelized counterpart relying on spectral multimode quantum optics. Section \ref{sec_pargen} is devoted to the investigation of parallelized generation of cat states and GKP codes. Finally, in Sec.~III we show that our protocol can be used with experimental readily accessible ressources such as single photon Fock states.

\section{Multimode iterative protocol}
\label{sec_prot}
\subsection{Framework: iterative protocols}
The protocol proposed in this paper takes its essence from iterative generation protocols, an elementary brick of which is represented in Fig. \ref{protelem}.\\
\begin{figure}[!h]
\begin{center}
\includegraphics[width=6cm]{Fig1}
\caption{Elementary brick of an iterative protocol.}
\label{protelem}
\end{center}
\end{figure}

Such operation consists in mixing two input `resource' states, $\ket{\psi_0}$ and $\ket{\psi_1}$, on a beamsplitter with transmittivity $\tau_0$, to herald the generation of an ouput state $\ket{\psi_{\rm out}}$. The heralding process can be performed in different ways~\cite{Etesse2014,Eaton2019}, and even if the results derived here can easily be extended to the photon counting case, we focus in the following on homodyne conditioning. In this scenario, one of the two output arms of the beamsplitter is measured by a homodyne detector (along a certain quadrature $x_\theta$). When the measurement outcome is close enough to an assigned value ($x_\theta=x_0$ on Fig.\ref{protelem}), the target output state is successfully heralded. The strength of such a protocol is that the average photon number of the output state can be higher than the one of the individual input states. Accordingly, it can be used to `breed' complex photonic states by cascading the protocol using $\ket{\psi_{\rm out}}$ as one of the input of a subsequent operation.\\
Figure \ref{protcasc} shows such a configuration, where states of increasing size are combined with resource states; tree-like configurations are also commonly investigated, thanks to the parallelization capabilities they offer \cite{Etesse14}.
\begin{figure}[!h]
\begin{center}
\includegraphics[width=\columnwidth]{Fig2}
\caption{A cascaded architecture for an iterative protocol.}
\label{protcasc}
\end{center}
\end{figure}
This idea is at the core of the `cat breeding' scheme, proposed in \cite{Lund2004}, where balanced beamsplitter ($\tau=1/\sqrt{2}$) and $x_i=0$ homodyne conditioning are used \cite{Takeoka2007,Brask2010}. In this case, mixing two SCSs allows to create a single SCS with an amplitude increased by a factor $\sqrt{2}$. Iterating this protocol $N$ times then allows to create a SCS with an amplitude amplified by a factor $2^{N/2}$. An alternative manner has also been proposed, where the input cat states were replaced by single photons, allowing the production of squeezed SCSs, robuster to losses \cite{Etesse2014,LeJeannic2018}. By modifying the quadrature along which the homodyne conditioning is performed ($p=0$ homodyne conditionings), GKP codes can be generated all-optically \cite{Vasconcelos2010,Etesse14}. More generally, all superpositions of up to $N$ photons can be generated by properly adjusting the beamsplitter reflectivities and homodyne conditionings \cite{Etesse2014}.\\
Proof-of principle experimental validations are currently limited to single stage implementations due to scalability issues \cite{Etesse2015,Sychev17}: in bulk optics, increasing the number of stages by increasing the number of input state sources indeed leads to intractable complexities. Temporal multiplexing represents a promising alternative that is also investigated~\cite{Bouillard2019}, but could ultimately lead to a decrease in the overall generation rate of the states. Integrated optics, on the other hand, offers an attractive solution to address spatial scalability issues, but inevitable losses associated with such platforms should be carefully taken into consideration \cite{Alibart2016}. \\
We propose in this article to rely instead on spectral multiplexing, which offers an increased generation efficiency as well as noticeable experimental simplifications.

\subsection{Multiplexed single-stage generation protocol}
The main idea of the proposal, sketched in Fig.\ref{protmult}, consists in the following key ingredients. To begin with, a `seed' state $\ket{\psi_0}$ initially populates a spectral mode $\chi_0(\omega)$ in which the quantum engineering process is to be performed. In output of the protocol, we will focus on the state $\rho_{\rm out}$ that populates this same mode $\chi_0(\omega)$. Subsequently, we select a collection of states $|\psi'_0\rangle$, $|\psi'_1\rangle$, ... $|\psi'_{N-2}\rangle$ respectively in modes $\varphi_0(\omega)$, $\varphi_1(\omega)$, ... $\varphi_{N-2}(\omega)$ that are used to `feed' the protocol, so as to grow the overall size of the output state $\rho_{\rm out}$.  As will be discussed with a concrete example, the spectral modes $\varphi_i(\omega)$ have to overlap $\chi_0(\omega)$ in a clever way, so that each $\ket{\psi'_i}$ contributes in an efficient manner to the generation protocol.
This collection of states is merged in a single spatial mode and sent to the first input port (port $a$) of a beamsplitter of transmittivity $\tau$ so as to interfere with the seed state $\ket{\psi_0}$ sent to the other input port of this beamsplitter (port $b$). Importantly, a perfect spatial mode matching of the two states populating modes $a$ and $b$ must be guaranteed. After interaction, one output port of the beamsplitter is measured by a homodyne detector along mode $\chi_0(\omega)$, to perform the desired conditioning $x_\theta=x_0$, for successfully heralding the generation of a state $\rho_{\rm out}$ on the other output port.
\begin{figure}[!h]
\begin{center}
\includegraphics[width=\columnwidth]{Fig3}
\caption{The multiplexed generation protocol (see text for description).}
\label{protmult}
\end{center}
\end{figure}

A major advantage of this method is that it allows mixing $N$ states in a single heralding step, therefore drastically increasing the success probability of the protocol and reducing the experimental complexity as compared with the standard spatial multiplexing strategy, as will be quantitatively shown later.\\
The multimode description that we use complies with realistic implementations of gaussian state sources, such as, for instance, multimode squeezed vacuum sources \cite{Wasilewski06}. It is now indeed well established that spontaneous parametric down-conversion (SPDC) through three-wave mixing process produces a collection of squeezed vacuum states in separate orthogonal spectral modes  $\{\varphi_k(\omega)\}_{k\in\mathbb N}$, each with a different squeezing parameter \cite{Branczyk2010,RomanRodriguez2021}.\\
More generally, each spectral mode $\varphi_k(\omega)$ can be associated with a bosonic anihilation operator $\hat{A}_k$ :
\begin{align}
\hat{A}_k=\int\varphi_k^*(\omega)\hat{a}(\omega)d\omega,
\end{align}
$\{a(\omega)\}$ being the bosonic operators of individual spectral components  at the SPDC output.
When two different sets of modes $\{\varphi_k(\omega)\}_{k\in\mathbb N}$ and $\{\chi_l(\omega)\}_{l\in\mathbb N}$ are involved, one can express their `likeness' by using the overlap formula
\begin{align}
c_{k,l}=\int\varphi_k^*(\omega)\chi_l(\omega)d\omega.
\end{align}
This formula allows switching from one mode basis to the other thanks to
\begin{align}
\varphi_k(\omega)=\sum_lc^*_{k,l}\chi_l(\omega).
\label{eq_baschang}
\end{align}
We show in the next paragraph that the proposed protocol can be used to generate highly-fidel SCSs and GKP codes with increased success probabilities compared to standard iterative protocols.

\section{Paralellized generation of complex photonic states}
\label{sec_pargen}
\subsection{Parallel cat `breeding'}
\label{subsec_catbreed}
\begin{figure*}
\begin{center}
\subfigure[]{
\includegraphics[width=0.49\columnwidth]{Fig4a.jpg}}
\subfigure[]{
\includegraphics[width=0.49\columnwidth]{Fig4b.jpg}}
\subfigure[]{
\includegraphics[width=0.49\columnwidth]{Fig4c.jpg}}
\subfigure[]{
\includegraphics[width=0.49\columnwidth]{Fig4d.jpg}}
\subfigure[$N = 2$]{
\includegraphics[width=0.49\columnwidth]{Fig4e.pdf}}
\subfigure[$N = 3$]{
\includegraphics[width=0.49\columnwidth]{Fig4f.pdf}}
\subfigure[$N = 6$]{
\includegraphics[width=0.49\columnwidth]{Fig4g.pdf}}
\subfigure[$N = 8$]{
\includegraphics[width=0.49\columnwidth]{Fig4h.pdf}}
\caption{Fidelity of the output state with a cat state of amplitude $\gamma$ as a function of the beamsplitter reflectivity $r$, for (a) $N=2$, (b) $N=3$, (c) $N=6$ and (d) $N=8$ input cat states of amplitude $\alpha=3$. In each plot, a star spots the highest fidelity point, and the corresponding parameters are given in the insets. The Wigner functions corresponding to this maximal fidelity point are represented in (e), (f), (g), (h) respectively.}
\label{cat_reflect}
\end{center}
\end{figure*}
The previously described scheme can be applied to the generation of large amplitude SCSs, inspired by cat `breeding' iterative protocols. To perform such a generation, we simply consider here $N\geqslant2$ identical SCS input states of amplitude $\alpha$ in the $N-1$ input modes $\varphi_k(\omega)$:
\begin{align}
\forall k,\in[\![0,N-2]\!],\ \ket{\psi_k'}=\ket{\psi_0}=\frac{1}{\mathcal{N}}\left(\ket{\alpha}+\ket{-\alpha}\right),
\end{align}
where $\mathcal{N}=\sqrt{2\left(1+e^{-2|\alpha|^2}\right)}$ is a normalization factor. Note that we restrict the discussion to even cat states, but similar reasoning is also valid with odd cat states, as treated in the Suppl. Mat.. Then we set that each mode $\varphi_k(\omega)$ is of equal contribution in the generation process, namely $\forall k\in[\![0,N-2]\!]$:
\begin{align}
c_{k,0}=\int\varphi_k^*(\omega)\chi_0(\omega)d\omega=\frac{1}{\sqrt{N-1}}.
\label{part_overlap}
\end{align}
In this way, each input `feeding' cat $\ket{\psi_k'}$ will `breed' the input seed cat $\ket{\psi_0}$ by the same amount.
Subsequently, homodyne conditioning is performed along the $x$ quadrature (in-phase with the SCS) around $x=0$, within a given acceptance window $[-\Delta x/2,\Delta x/2]$. An analytic formula for the output Wigner function (see Suppl. Mat.) can then be derived, allowing us to predict the performances of the protocol as a function of key parameters such as $r$, $N$, $\alpha$ and the heralding width $\Delta x$. Many different metrics can be used to quantify the distance between an output state and a target state. Here, the fidelity of $\rho_{\rm out}$ with the nearest ideal SCS is used (amplitude $\gamma$).
\subsubsection*{Beamsplitter reflectivity}
The aforementioned fidelity of the output state is plotted in Figs. \ref{cat_reflect}(a) to (d) for $N=$2, 3, 6 and 8  input SCSs of amplitude $\alpha=3$, as a function of the beamsplitter reflectivity $\tau$ and of the output state amplitude $\gamma$. These plots clearly reveal that the output amplitude of the cat state in mode $\chi_0(\omega)$ that corresponds to optimal fidelities with $\rho_{\rm out}$ grows with the number of input states. This successfully proves the breeding of the input SCSs. Figure~\ref{cat_reflect} also shows that the reflectivity should be optimized in a non-trivial way, according to the number of input breeding states, in order to maximize the output state fidelity with high-amplitude SCSs.

To gain a qualitative intuition on the quality of the multimode breeding operation, the Wigner functions of the corresponding output states are plotted in Figs.\ref{cat_reflect} (e) to (h), showing their strong non-classical signature, as the Wigner functions contain oscillations with negative values. To be more quantitative on this result, the numerical values of fidelity, optimal reflectivity and output SCS amplitudes are also written on Figs. \ref{cat_reflect} (a) to (d), revealing the very high fidelities of the output states with high-amplitude SCSs. In addition, it confirms that the amplitude of the output SCS increases with $N$, with a scaling of the order of $\gamma~\sim~\sqrt{N}\alpha$.
\subsubsection*{Amplitude of the input states}
The protocol is also sensitive to the amplitude of the input cat states. This is illustrated in Fig.~\ref{dep_alpha}, where the output state fidelity with the nearest cat state is represented as a function of the input SCS amplitude $\alpha$, for different numbers of input states $N$. Two important trends shall be noticed: (i) for a given $\alpha$, increasing $N$ decreases the overall output state fidelity. However, (ii) increasing $\alpha$ for a given $N$ leads to an increased global fidelity. For instance, to create states with 99\% fidelity with $N=2$, an input amplitude $\alpha=1.13$ is sufficient, whereas $\alpha=2.49$ is required to reach the same fidelity with $N=4$ input states.
 In the case of even cat states that are considered here, there is also an increase of the fidelity for vanishing $\alpha$ simply because vacuum is a perfect even cat state with zero amplitude, a case of no interest for us.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.9\columnwidth]{Fig5}
\caption{Fidelity of the output state with the nearest SCS as a function of $\alpha$, for different $N$. The cases represented in fig.\ref{cat_reflect} correspond to the vertical blue line ($\alpha=3$).}
\label{dep_alpha}
\end{center}
\end{figure}

\subsubsection*{Comparison with the iterative cat breeding protocol}
Interestingly, we can draw a direct parallel with the iterative cat breeding protocol, where $N-1$ heralding steps are performed on $N$ input resources. The main difference, however, is that the multimode protocol that we propose here is inevitably accompanied with information loss due to partial overlap of the spectral modes $\varphi_k(\omega)$ with the output mode of interest $\chi_0(\omega)$. Equation (\ref{part_overlap}) directly shows that each mode $\varphi_i(\omega)$ overlaps only partially with $\chi_0(\omega)$, therefore all of the components of the input states $\ket{\psi'_i}$ along $\chi_{l>0}(\omega)$ are purely lost during the process, where $\{\chi_l(\omega)\}_{l\in\mathbb N}$ constitutes a complete spectral basis. Therefore, the output state $\rho_{\rm out}$ cannot be a pure state as long as more than two modes are populated in total.\\
To give an intuitive model, we propose a simple parallel between spatial and spectral modes: the operation depticted in Fig.~\ref{protmult} by the merging arrows from modes $\varphi_k(\omega)$ to $\chi_0(\omega)$ can actually be simply modeled as a $N-1$ input / $N-1$ output beamsplitter transformation, where $N-2$ outputs are lost in the environment. This transformation is exactly the one depicted in Fig.~\ref{protcasc} for input states $\ket{\psi_0}$ to $\ket{\psi_{N-2}}$, with $\tau_i=1/\sqrt{i+2}$ so as to obtain $c_{0,k}=1/\sqrt{N-1}$. In other words, our protocol is equivalent to the iterative breeding scheme of Fig.~\ref{protcasc}, where only the last heralding step $x=x_{N-1}$ is performed and all the other ones are simply ignored. The direct consequence for this is that the overall success probability of the protocol is increased by several orders of magnitude as compared to the iterative scheme. This is illustrated in Fig.~\ref{psucc_compare}, where the success probability of our scheme for input SCS amplitudes of $\alpha=3$ and $N=4$ is plotted, in comparison with the iterative scheme of Fig.\ref{protcasc} in the case where all the heralding widths are identical. The cost obviously is the overall fidelity, which is degraded in the multimode case, but nevertheless remains very high, as we have seen in Fig.~\ref{cat_reflect}. The multimode protocol therefore takes advantage of the robustness of the breeding protocol for imperfect heralding at sufficiently high SCS amplitudes.
\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.9\columnwidth]{Fig6}
\caption{Comparison between the success probability of the multimode scheme and the iterative scheme for $N=4$, as a function of the heralding window $\Delta x$ normalized to the shot noise standard deviation $\sigma_0$.}
\label{psucc_compare}
\end{center}
\end{figure}
\subsection{GKP states preparation}
In addition to the generation of SCSs, our protocol can also advantageously be applied to for the generation of more complex states such as the GKP codes \cite{Gottesman2001}. These states were first introduced as a means to encode information in an oscillator, with an error-correcting capability. They were recently generated and used in trapped ions~\cite{Fluhmann2019} and superconducting circuits~\cite{CampagneIbarcq2020}, but still remain out of reach of all-optical protocols. Here, we can use the same protocol as the one depicted in Sec. \ref{subsec_catbreed} ($N-1$ input cat states in modes $\varphi_k$ and one cat state in mode $\chi_0$, all of amplitude $\alpha$), with the only difference that homodyne conditioning is performed along $x_{\pi/2}:=p=0$. Such conditioning has already been investigated in the context of all-optical generation of GKP with iterative protocol~\cite{Vasconcelos2010,Etesse14}. We show here that it can be transposed to our multiplexed approach. As a reminder, the wavefunctions of approximate GKP codes $\ket{\tilde{0}}$ and $\ket{\tilde{1}}$ take the form
\begin{subequations}
\begin{align}
\langle x \vert \tilde{0} \rangle& = \dfrac{1}{\mathcal{N}_0} \sum_{k \in \mathbb{Z}} e^{-\frac{s_1^2}{2} (ka)^2} e^{-\frac{1}{2 s_2^2}(x - k a)^2}\\
\langle x \vert \tilde{1} \rangle &= \dfrac{1}{\mathcal{N}_0} \sum_{k \in \mathbb{Z}} e^{-\frac{s_1^2}{2} (k + 1 / 2)^2 a^2} e^{-\frac{1}{2 s_2^2}\left(x - (k + 1 / 2) a \right)^2},
\end{align}
\end{subequations}
where $1/s_1$ gives the overall envelope size of the state, $s_2$ is the standard deviation of each tooth in the distribution, and $a$ is the distance between consecutive teeth along the $x$ quadrature.
\begin{figure}[!h]
\begin{center}
\subfigure[]{
\includegraphics[width=0.8\columnwidth]{Fig7a.png}}
\subfigure[]{
\includegraphics[width=0.8\columnwidth]{Fig7b.pdf}}
\caption{Performance of the multiplexed protocol for the generation of approximate GKP codes. (a) Wigner function of the output state with $N=3$ input SCSs of amplitude $\alpha=3$. (b) GKP state amplitude when the amount of input cat states $N$ is augmented.}
\label{GKP_perfs}
\end{center}
\end{figure}
Figure \ref{GKP_perfs}(a) shows the shape of the Wigner function of the output state for $N=3$, revealing the grid structure specific to GKP codes. Such likeness is, here again, quantified by the fidelity with the nearest GKP code, that we plot in Fig. \ref{GKP_perfs}(b). This figure shows that augmenting the amount of input resources indeed increases the size of the output state along the opposite quadrature ($1/a$), while requiring an optimization of $r$. Similarly to the SCS case, it also shows that the fidelity decreases with $N$.

\section{Proposal with experimentally accessible resources}
\begin{figure*}
\begin{center}
\subfigure[]{
\includegraphics[width=0.45\columnwidth]{Fig8a.png}}
\subfigure[]{
\includegraphics[width=0.45\columnwidth]{Fig8b.png}}
\subfigure[]{
\includegraphics[width=0.45\columnwidth]{Fig8c.png}}
\subfigure[]{
\includegraphics[width=0.45\columnwidth]{Fig8d.png}}

\caption{Wigner functions of the output state of the protocol fed by single photon Fock states, with (a) $N=2$, (b) $N=3$, (c) $N=4$ and (d) $N=5$.}
\label{sing_phot}
\end{center}
\end{figure*}
An important aspect of our strategy is that it can directly be implemented with already accessible experimental resources. Even if high-amplitude SCSs are not trivial to generate experimentally, good approximation of them at low amplitudes (also called Schr\"odinger `kitten' states) are commonly produced states. Moreover, the breeding principle has been experimentally demonstrated with kittens, showing the validity of the approach \cite{Sychev17}. In addition, single photon Fock states also constitute a valid alternative to them, as theoretically and experimentally demonstrated \cite{Etesse14,Etesse2015}, to generate squeezed SCS-like output states with high fidelity. It is worth noting that squeezing the output states constitutes a high advantage, as it can increase their robustness to losses~\cite{Ourjoumtsev2006,LeJeannic2018}. Accordingly, the results presented here will focus on input single photon Fock states. Very similar results are expected with kitten states.\\
We recall that our protocol requires the use of a multiplexed input state source. This can simply be achieved by taking advantage of the spectral tailoring possibility that processes like spontaneous parametric down-conversion offer~\cite{Branczyk2010}. Indeed, depending on the pumping and phase-matching conditions, two-mode squeezing can be generated along a large number of spectral modes, offering the ideal platform for multiplexed generation~\cite{RomanRodriguez2021}. Then, single photon Fock states are simply heralded by performing a projective measurement over one of the two mode squeezing subsystems. This method, although very commonly used for the generation of high-quality single photon Fock states for a couple of decades now, has recently been extended experimentally to high-dimension mode-selective quantum state tailoring~\cite{Ra2020}. This allows to tweak the supermode in which the single photon detection occurs, therefore heralding a single photon Fock state in a given supermode. Here we can take advantage of such capabilities, by using multiplexed photon sources to feed the protocol and generate complex photonic states in a single shot.\\
The output Wigner functions of the protocol fed by an increasing number $N$ of single photon Fock states are shown in Fig.~\ref{sing_phot}. They clearly reveal that the output state contains more and more oscillations at its center with increasing $N$. To be more quantitative, Tab. \ref{table_perf_SP} gathers the parameters of the closest SCS, where $\alpha$ is the amplitude, $s$ its compression parameter and $\mathcal{F}$ the fidelity between the two states. This table reveals that the amplitude of the state increases with $N$ for low values, but degrades after $N=4$. Further simulations show that this is true for larger $N$ as well. Such an evolution can qualitatively be explained by the phase invariance of the two states impinging the beamsplitter: as the only rotation symmetry breaking operation is the homodyne measurement, the fidelity of the created state with a highly non-phase invariant state decreases with its size.
\begin{table}[h!]
\begin{tabular}{|m{1cm}||m{1cm}|m{1cm}|m{1cm}|m{1cm}|}
\hline
&$N=2$&$N=3$&$N=4$&$N=5$\\
\hline
\hline
$\alpha$&1.63&2.03&2.30&2.12\\
\hline
$s$&0.66&0.65&0.66&0.74\\
\hline
$\mathcal{F}$&0.99&0.90&0.66&0.49\\
\hline
\end{tabular}
\caption{Fidelity of the output state with a squeezed SCS of amplitude $\alpha$ and squeezing $s$.}
\label{table_perf_SP}
\end{table}
It should also be pointed that the optimal generation condition (highest amplitude and fidelity) is fulfilled with the use of a low-reflectivity beamsplitter (see Suppl. Mat. for numerical simulations).\\
This analysis reveals that few-mode single photon states are relevant resources for our protocol in order to produce high-quality SCSs, essential brick in continuous-variable quantum information processing.

\section*{conclusion} We have proposed and discussed an innovative protocol for the generation of non-Gaussian states of light. Our scheme is inspired by iterative protocols but takes advantage of multiplexing capabilities of photonics resources. It relies on the use of multiple states in a single spatial mode but multiplexed in an other degree of freedom, mixed with a resource state whose mode partially overlaps all input states. By performing a homodyne conditioning after mixing the states at a beamsplitter, complex output states can be produced. The major strength of our proposal is that no cascade manipulation steps are required for the state breeding, therefore strongly simplifying the experimental implementation in a one-shot scheme, and allowing to access high-success rates. We have shown that frequency multiplexing can be used for the generation of high-quality and high-amplitude Schr\"odinger cat and GKP states, by using as input states small-amplitude Schr\"odinger cat states, or alternatively, readily accessible experimental resources such as multiplexed single photon Fock states. Note that current state-of-the-art non-linear sources offer the capability for generating such states by providing the native generation of frequency multiplexed single- and two-mode squeezing. Due to its high efficiency and simplified implementation as well as to the possibility of implementing it with already available experimental resources, we believe that our approach to non-Gaussian state preparation represents a major breakthrough in the domain of non-classical photonic state generation, and will open a wide range of applications.
\section*{acknowledgments}
This work has been conducted within the framework of the project HyLight (No. ANR-17-CE30-0006-01) granted by the Agence Nationale de la Recherche (ANR).