\section{Introduction}
\label{sec:introduction}

Interstellar magnetic fields (\emph{B}-fields) are thought to play an important role in the formation of molecular clouds.
As material from the interstellar medium (ISM) flows along field lines onto the nascent clouds, the resulting filamentary structures are expected to extend perpendicular to the interstellar \emph{B}-field (\citealp{2019FrASS...6....5H} for a review).
However, this flow of matter can also influence the morphology of the surrounding \emph{B}-field.
Therefore, the \emph{B}-field structures we observe around molecular clouds are imprinted with their history of an accumulation from the ISM, thus allowing us to trace the processes that led to their formation \citep[e.g.,][]{2018MNRAS.480.2939G}.

The plane-of-sky (POS) component of the \emph{B}-field can be traced by polarimetric observations of optical and near-infrared radiation from stars located behind molecular clouds, as well as thermal continuum emission from interstellar dust particles in those same clouds \citep{1993prpl.conf..279H, 2007JQSRT.106..225L, 2008MNRAS.388..117H, 2009ApJS..182..143M, 2012ARA&A..50...29C}.
Aspherical dust particles irradiated by starlight are charged up by the photoelectric effect, as well as spun up as a result of radiative torques (RATs; \citealp[e.g.,][]{1996ApJ...470..551D,1997ApJ...480..633D,2019ApJ...883..122L}).
These spinning particles are aligned with their rotation axes (i.e., their minor axes) parallel to the \emph{B}-field orientation.
This alignment of the dust particles results in preferential absorption and scattering of the background starlight, which makes the observed starlight polarized in the direction parallel to the POS component of the \emph{B}-field.
In the case of thermal dust emission in sub-mm wavelengths, the preferential emission from the aligned dust particles causes the emission polarized in the direction perpendicular to the POS component of the \emph{B}-field \citep{1966ApJ...144..318S, 1988QJRAS..29..327H, 2015ARA&A..53..501A}.

The Perseus molecular cloud, whose distance is about 300 pc from the sun \citep{2018ApJ...865...73O, 2018ApJ...869...83Z, 2021A&A...645A..55P}, is one of the most active star-forming molecular clouds in the solar vicinity \citep{2008hsf1.book..308B}.
This cloud is associated with a single H\small{I} shell along with the Taurus, Auriga, and California molecular clouds.
It has been suggested that this H\small{I} shell may have been formed by the interstellar bubble around the Perseus OB2 association \citep{2008hsf1.book..308B, 2013ApJ...765..107L, 2019A&A...623A..16S}.

\citet{2020ApJ...899...28D} revealed that the active star-forming region NGC 1333 in the Perseus molecular cloud shows a complex \emph{B}-field structure at spatial scales $< 1$ pc.
In their analysis, the small scale variation of the \emph{B}-field is well explained if the \emph{B}-field is generally perpendicular to the local dense ISM filaments.
Fluctuations in the \emph{B}-field structure at small spatial scales in the Perseus molecular cloud are also found by optical polarimetry (\citealp{1990ApJ...359..363G}; Figure \ref{fig:Goodman}).
The observed polarization angles indicate that POS orientations of the \emph{B}-field ($\theta_\mathrm{star}$) show two distinct populations.
One population has a peak at $\theta_\mathrm{star} \sim -40 ^\circ$ and the other population has a broader peak in position angles at $\theta_\mathrm{star} \sim +70 ^\circ$.
The two populations of vectors show no spatial segregation across the Perseus cloud complex (\citealp{1990ApJ...359..363G}; see Figure \ref{fig:Goodman}).

It is not yet known whether these bi-modal $\theta_\mathrm{star}$ values represent the local variation of \emph{B}-field orientations, differences in the characteristics of the dust particles, or whether that is caused by the superposition of multiple ISM components along the line of sight \citep[LOS;][]{1990ApJ...359..363G, 2002ApJ...574..822M, 2006ApJ...643..932R, 2019ApJ...871L..15G}.
In this paper, we combine optical polarimetry data with Gaia measurements of stellar distances, to determine the multi-layer distribution of the \emph{B}-field in the direction of the Perseus molecular cloud.
With this analysis, we aim to reveal the cause of the bi-modal \emph{B}-field structure observed in the Perseus molecular cloud region, which we find due to a contribution of the foreground Taurus molecular cloud.

This paper is organized as follows.
In Section \ref{sec:data}, we describe the data used for our analysis.
In Section \ref{sec:3D}, we analyze the stellar distances of optical polarimetry and identify contributions from the Perseus and the foreground Taurus molecular clouds.
In Section \ref{sec:discussion}, we discuss the relationship between the H\small{I} shell, which is thought to be associated with the Perseus-Taurus molecular cloud, and the \emph{B}-field distribution we have identified.
We also discuss the alignment between $\theta_\mathrm{star}$ and the Planck-observed \emph{B}-field orientation ($\theta_\mathrm{Planck}$), which gives information on small-scale \emph{B}-field structure below Planck's beam size.
In Section \ref{sec:conclusion}, we summarize the results.

\section{Data}
\label{sec:data}

\subsection{Optical and Near-Infrared Stellar Polarimetry}

We use optical polarimetry of 88 sources in the direction of the Perseus molecular cloud observed by \citet[][Figure \ref{fig:Goodman}]{1990ApJ...359..363G}.
The passband of the observations were centered at 762.5 nm with a bandwidth of 245 nm.

We also use near-infrared (NIR) polarimetry data taken toward NGC 1333 in the Perseus molecular cloud, as shown in the inset of Figure \ref{fig:Goodman}.
We use \emph{K}-band polarimetry by \citet[][14 sources]{1988MNRAS.231..445T} and \emph{R}- and \emph{J}-band polarimetry by \citet[][33 sources]{2011AJ....142...33A}.
The near-infrared $\theta_\mathrm{star}$ are mainly distributed from $\theta_\mathrm{star} \sim -90^\circ$ to $\theta_\mathrm{star} \sim -40^\circ$.

\subsection{Planck Sub-mm Polarimetry}
\label{sec:Planck}

We estimate the \emph{B}-field orientation measured in sub-mm dust emission by using Planck data at 353 GHz \citep{2020A&A...641A..12P}.
We use \verb|HFI_SkyMap_353-psb_2048_R3.01_full.fits| taken from the Planck Legacy Archive, https://pla.esac.esa.int/.
In our analysis, we set the spatial resolution of the Planck polarimetric data as a $10'$ FWHM Gaussian to achieve good S/Ns.

\subsection{Gaia DR2 Photometry and Trigonometric Distances}
\label{sec:DR2}

We estimate the stellar distances by using Gaia astrometry data (DR2; \citealp{2016A&A...595A...1G, 2018A&A...616A...1G}), and use SIMBAD \citep{2000A&AS..143....9W} positions to cross-match the stars in the Gaia catalog.
We set a search radius of $5\arcsec$ and take the star at the closest position for \emph{R}-band and \emph{J}-band data by \citet{2011AJ....142...33A}.
The identified stars show $M_G = 14$ -- 20 mag, and have good correspondence with the \emph{J}-band magnitude tabulated by \citet{2011AJ....142...33A}.

For relatively old datasets by \citet{1990ApJ...359..363G} and \citet{1988MNRAS.231..445T}, we set a relatively large search radius of $30\arcsec$ and take the brightest star in the searched region, which gives a good cross-match of the catalogued stars as follows.
We find that the stars observed by \citet{1990ApJ...359..363G} show \emph{G}-band magnitude in the Gaia catalog between $M_G = 7$ -- 15 mag.
We exclude three stars that show $M_G = 17$ -- 21 mag from the following analysis for their possible misidentifications.
The stars observed by \citet{1988MNRAS.231..445T} show $M_G = 10$ -- 18 mag. We exclude one star that shows $M_G = 20$ mag from the following analysis because of its non-reliable distance estimation (negative parallax).

The Gaia parallax in DR2 is known to have a systematic bias of -0.03 mas (see \citealp{2015PASP..127..994B}; \citealp{2018AJ....156...58B}; \citealp{2018A&A...616A...2L}; \citealp{2018A&A...616A..17A}; \citealp{2019A&A...621A..48L}; \citealp{2020Ap&SS.365..112M}).
This bias is negligible in our distance evaluation, as the stellar distances in our analysis are less than 1 kpc and thus their parallaxes $\varpi > 1$ mas.
Thus, we did not correct the bias, and estimate the distances of each star by $d = 1000/\varpi$ (pc).

The Renormalised Unit Weight Error (RUWE) is a parameter that is expected to be around 1.0 for sources where the single-star model provides a good fit to the astrometric observations.
A value significantly greater than 1.0 (say, $>1.4$) could indicate that the source is non-single or otherwise problematic for the astrometric solution
(see the Gaia data release documentation\footnote{https://gea.esac.esa.int/archive/documentation/GDR2/}).
We estimate the RUWE for each source following the formulation described in the document ``Re-normalising the astrometric chi-square in Gaia DR2”\footnote{https://www.cosmos.esa.int/web/gaia/public-dpac-documents} and use the data only if whose $\mathrm{RUWE} \leq 1.4$.

As a result, we estimate the distances of 111 sources, including 70 sources from \citet{1990ApJ...359..363G}, 20 sources (\emph{R}-band) and 15 sources (\emph{J}-band) from \citet{2011AJ....142...33A}, and 6 sources from \citet{1988MNRAS.231..445T}.
We summarize the identified data in Table \ref{tab:stellar_data} in Appendix \ref{sec:used_data}.

\section{Results}
\label{sec:3D}

\subsection{Distance to the Clouds that Produce Polarization}
\label{sec:breakpoints}

\begin{figure*}[tp]
\centering
\includegraphics[width=\linewidth]{Perseus_stardist.png}
\caption{Stellar polarimetry data toward the Perseus molecular cloud (black line segments), overlaid on the Planck-observed \emph{B}-field orientation \citep[gray line segments;][]{2020A&A...641A..12P}.
Black line segments in the main panel are the optical polarimetry data by \citet{1990ApJ...359..363G}, and those in the inset are near-infrared polarimetry data by \citet{1988MNRAS.231..445T} and \citet{2011AJ....142...33A} toward NGC 1333.
The length of the black line segments is proportional to the polarization fraction ($P_\mathrm{star}$).
Reference scales of $P_\mathrm{star}$ are shown in the lower right corner of both the main panel and the inset.
Filled circles indicate the position of the stars and are color-coded to indicate stellar distance ranges observed with Gaia.
Group 1 is for distances $d > 300$ pc, Group 2 is for $150 < d < 300$ pc, and Group 3 is for $d < 150$ pc.
N/A is for the data with no reliable distance estimation available.
See Section \ref{sec:breakpoints} for the discussion and the definition of these distance groups.
For the Planck data (gray line segments), the line segments' length has been normalized to show only the orientation of the \emph{B}-field.
The spatial resolution of the Planck data is $10'$.
The background color scale is the hydrogen column density ($N_\mathrm{H}$).
We convert the dust opacity at 353 GHz ($\tau_{353}$), estimated by Herschel and Planck observational data \citep{2016A&A...587A.106Z}, to the hydrogen column density  by $N_\mathrm{H} = 1.6\times 10^{26}~\tau_{353}$ (cm$^{-2}$; \citealp{2014A&A...571A..11P,2015A&A...576A.104P}).
A reference scale of 2 pc, in which we assume the distance to the field is 300 pc, is shown in the lower left corner of the main panel.
}
\label{fig:Goodman}
\end{figure*}

Figure \ref{fig:Goodman} compares the spatial distribution of the polarimetry data in our analysis with the inferred \emph{B}-field morphology from Planck.
The polarization position angles, $\theta_\mathrm{star}$, show a bi-model distribution, with concentrations at $\theta_\mathrm{star} \sim +70^\circ$ and $\theta_\mathrm{star} \sim -40^\circ$ (measured from North to East).
\citet{1990ApJ...359..363G} found that the two populations in $\theta_\mathrm{star}$ also have distinct polarization fractions ($P_\mathrm{star}$).
We display the relationship between the observed $\theta_\mathrm{star}$ and $P_\mathrm{star}$ in Figure \ref{fig:GoodmanAngP}.
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{Goodman_ANG_P}
\caption{
A correlation between $\theta_\mathrm{star}$ and $P_\mathrm{star}$ of the optical polarimetry by \citet{1990ApJ...359..363G}.
Histgrams of the two parameters are also shown on each axis.
Color codes are the stellar distance ranges same as Figure \ref{fig:Goodman}.
Yellow data points are the data with no reliable distance estimation available.
}
\label{fig:GoodmanAngP}
\end{figure}
The population of $\theta_\mathrm{star} \sim -40^\circ$ show relatively larger $P_\mathrm{star} \gtrsim 1$\%, while the other population that has $\theta_\mathrm{star} \sim +70^\circ$ show relatively smaller $P_\mathrm{star} \lesssim 1$\%.

In Figure \ref{fig:AngPbreak}, we display the $\theta_\mathrm{star}$ and $P_\mathrm{star}$ dependences as a function of the estimated stellar distances.
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{ANG_P_breaks}
\caption{The relationships between the stellar distances $d$ and $\theta_\mathrm{star}$ and $P_\mathrm{star}$.
Note that the horizontal axis is a linear scale below 300 pc, but a logarithmic scale above 300 pc.
Purple points are the optical polarimetry by \citet{1990ApJ...359..363G}.
Blue points are the \emph{R}-band polarimetry and green points are the \emph{J}-band polarimetry, respectively, by \citep{2011AJ....142...33A}.
Yellow points are \emph{K}-band polarimetry by \citet{1988MNRAS.231..445T}.
Three stellar groups based on their distances (Group 1, Group 2, and Group 3: see Section \ref{sec:breakpoints} and Table \ref{tab:ANG_P_Group}) are indicated.
Vertical dashed lines indicate the breakpoints in the distribution of the data points, and shaded area between dotted lines indicate their 95\% confidence intervals.
See Section \ref{sec:breakpoints} for the detailed description of the breakpoints estimation.
The horizontal dashed line in the top panel shows the orientation parallel to the galactic plane.
}
\label{fig:AngPbreak}
\end{figure}
As seen in the figure, there is a clear jump in both distributions at a distance of about 300 pc, which is the distance to the Perseus molecular cloud.
In addition to that, another jump is noticeable at a distance of about 150 pc.
Hereafter, we call the polarimetry data whose stellar distances $d > 300$ pc as Group 1, those with $150 < d < 300$ pc as Group 2, and those with $d < 150$ pc as Group 3, respectively.
We summarize $\theta_\mathrm{star}$ and $P_\mathrm{star}$ values of each group in Table \ref{tab:ANG_P_Group}.
\begin{table}[t]
\begin{center}
\caption{$\theta_\mathrm{star}$ and $P_\mathrm{star}$ Values of Each Distance Group obtained from Optical Polarimetry Data \citep{1990ApJ...359..363G}}
\label{tab:ANG_P_Group}
\begin{tabular}{lcccc}
\hline
\hline
& Distance & $N_\mathrm{star}$$^\mathrm{a}$ & $\theta_\mathrm{star}$$^\mathrm{b}$ & $P_\mathrm{star}$\\
& (pc) & & (deg) & (\%)\\
\hline
Group 1 & $> 300$ & 39 & $-37.6 \pm 35.2$ & $2.0 \pm 1.7$ \\
Group 2 & 150 -- 300 & 25 & $+66.8 \pm 19.1$ & $0.8 \pm 0.6$ \\
Group 3 & $< 150$ & 6 & $-41.8 \pm 32.0$ & $0.1 \pm 0.1$ \\
\hline
\multicolumn{5}{p{0.44\textwidth}}{$^\mathrm{a}$The number of stars.}\\
\multicolumn{5}{p{0.44\textwidth}}{$^\mathrm{b}$$\theta$'s mean and standard deviation values are circular means and circular standard deviations, which take into account the $\theta$'s $180^\circ$ degeneracy, throughout this paper.
The definitions of the circular mean and the circular standard deviation are given by \citet{2020ApJ...899...28D}.}
\end{tabular}
\end{center}
\end{table}

Figure \ref{fig:AngPbreak} shows that $\theta_\mathrm{star}$ and $P_\mathrm{star}$ are both consistent within their own distance group.
That is, the stars in Group 1 and the stars in Group 2 are tracing the polarization from ISM clouds located at distances of 300 pc and 150 pc, respectively.
In particular, the consistency of $P_\mathrm{star}$ indicates that the ISM between and behind the two clouds does not significantly contribute to stellar polarization.
Group 3 can be thought as foreground stars with little interstellar extinction.
We note that these stars have very high uncertainties in their polarization measurements due to having very low polarization fractions.

To quantitatively estimate the distance of the two ISM clouds that cause polarization, we perform a breakpoint analysis on $\theta_\mathrm{star}$ and $P_\mathrm{star}$ distributions as a function of $d$ shown in Figure \ref{fig:AngPbreak}.
We assume that $\theta_\mathrm{star}$ and $P_\mathrm{star}$ are constant as a function of $d$, which corresponds to the assumption that the observed polarization is caused by 2D sheet(s) of ISM at specific distance(s).
In addition, we assume the $\theta_\mathrm{star}$ and $P_\mathrm{star}$ distributions have a certain number of step-wise changes (i.e., breakpoints), which correspond to the positions of the 2D sheets.
We perform least-squares fits to the data and make most likelihood estimations (MLE) of the positions of breakpoints.
We then repeat the fit with different number of breakpoints, and compare the goodness-of-fit values based on the Bayesian information criterion, to get the most likely numbers of breakpoints and their positions.

We perform the breakpoint analysis for each distance dependence of $\theta_\mathrm{star}$ and $P_\mathrm{star}$ shown in Figure \ref{fig:AngPbreak} by using the R library `strucchange' \citep{strucchange,breakpoints}.
For $\theta_\mathrm{star}$ distribution analysis, we combine optical and NIR polarimetry data as shown in Figure \ref{fig:AngPbreak} and use the shortest waveband data if multiple waveband data are available for a single star.
On the other hand, we use only optical polarimetry data to analyze $P_\mathrm{star}$ distribution because different wavelengths give different polarization fractions.
We used $P_\mathrm{star}$'s logarithm values in our analysis to detect breakpoints of widely different $P_\mathrm{star}$ values (see Figure \ref{fig:AngPbreak}) with comparable sensitivity to each other.

The estimated distances of the breakpoints are shown in Figure \ref{fig:AngPbreak} and Table \ref{tab:ANG_P_breaks}.
The breakpoint analysis shows that there are two breakpoints in both $\theta_\mathrm{star}$ and $P_\mathrm{star}$ distributions.
The estimated distances are $150^{+4}_{-34}$ pc and $303^{+12}_{-7}$ pc for $\theta_\mathrm{star}$, and $150^{+4}_{-136}$ pc and $295^{+22}_{-67}$ pc for $P_\mathrm{star}$, respectively.
Here we show the MLE positions and their 95\% confidence intervals.
Although the breakpoints in $P_\mathrm{star}$ have larger errors comparing to that in $\theta_\mathrm{star}$ due to smaller jumps in $P_\mathrm{star}$ values, the estimated distances are fully compatible between $\theta_\mathrm{star}$ and $P_\mathrm{star}$.
\begin{table}[t]
\begin{center}
\caption{Estimated Cloud Distances as $\theta_\mathrm{star}$ and $P_\mathrm{star}$ Breakpoints}
\label{tab:ANG_P_breaks}
\begin{tabular}{lll}
\hline
\hline
& \multicolumn{2}{c}{Cloud Distance}\\
& \multicolumn{2}{c}{(pc)}\\
\hline
$\theta_\mathrm{star}$ & $150^{+4}_{-34}$ & $303^{+12}_{-7}$ \\
$P_\mathrm{star}$ & $150^{+4}_{-136}$ & $295^{+22}_{-67}$ \\
\hline
\end{tabular}
\end{center}
\end{table}

The distance of the breakpoint at 300 pc is consistent with the distance to the Perseus molecular cloud \citep{2018ApJ...865...73O, 2018ApJ...869...83Z,2021A&A...645A..55P}.
Thus, we conclude that the polarization of Group 1 traces the \emph{B}-field of the Perseus molecular cloud.
The distance of another breakpoint at 150 pc is thought to indicate the distance to the foreground ISM cloud that causes the polarization of Group 2.
This 150 pc foreground ISM also contribute to depolarize Group 1.
We will discuss this depolarization effect in Section \ref{sec:qu_separation}.
Group 3 are the foreground stars that show no clear indication of interstellar extinction.

\subsection{Estimation of Cloud Distances based on Photometry and Gaia Distances}
\label{sec:Gaia_Phot}

In the previous section (Section \ref{sec:breakpoints}), we estimated the distances of the breakpoints found in polarimetric data.
The accuracy of estimation is limited to a few tens of parsecs because of small number of data points.
In this section, we perform a breakpoint analysis using all the available stellar photometric data in the Perseus molecular cloud direction to make a more accurate estimation of the distances to the breakpoints, including the foreground ISM.
We then compare our result with other existing estimates of cloud distances in this direction to deomonstrate the reliability of our method for measuring dust clouds' distance.

We use \emph{G}-band extinction ($A_G$) values catalogued in Gaia DR2 with $d \leq 1$ kpc and $\mathrm{RUWE} \leq 1.4$, and fit these values as a function of distances obtained from Gaia parallax.
The spatial distribution of the stellar data used in this analysis is shown in Figure \ref{fig:PerseusGaiadist}.
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{Perseus_Gaia_dist.png}
\caption{
Distribution of Gaia photometric data in the direction of the Perseus molecular cloud.
We plot the stars whose distances $d \leq 1$ kpc and $\mathrm{RUWE} \leq 1.4$ (see Section \ref{sec:DR2} for the definition of RUWE).
The background color scale and contours are the hydrogen column density estimated by Herschel and Planck observational data \citep[][see Figure \ref{fig:Goodman}]{2016A&A...587A.106Z}.
Contour levels are $N_\mathrm{H} = 0.5,~1,~\mathrm{and}~5 \times 10^{22} ~\left( \mathrm{cm}^{-2} \right)$.
Red circles (diameter: $1^\circ$) are the regions where we estimate the breakpoint distances and compare them with the CO molecular cloud distances estimated by \citet{2018ApJ...869...83Z}.
Names of the CO molecular clouds are shown above the figure.
See Section \ref{sec:Gaia_Phot} and Figure \ref{fig:COMPLETE_distance} for the details of the comparison.}
\label{fig:PerseusGaiadist}
\end{figure}
We note the paucity of stars in the direction of the dense molecular cloud.
This is due to the larger extinction in the visible \emph{G}-band compared with the NIR bands.

\citet{2018ApJ...869...83Z} combined Gaia distance and NIR photometry to obtain the distance of each velocity component of the $^{12}\mathrm{CO}~(1-0)$ line emission for each cloud core in the Perseus molecular cloud.
We perform breakpoint analyses for similar spatial directions and compare the result with their estimation.
Since the exact spatial region analyzed in \citet{2018ApJ...869...83Z} is not indicated in the literature, we estimate the breakpoints at $1^\circ$ diameter regions centered at each molecular cloud core.
The regions are shown as red circles in Figure \ref{fig:PerseusGaiadist}.

\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{COMPLETE_distance}
\caption{
Estimated distances of major star-forming regions across the Perseus molecular cloud.
Results of our breakpoint analysis (`this work' in the figure) are compared with the CO distances by \citet[][also see Table \ref{tab:COMPLETE_distance}]{2018ApJ...869...83Z}.
Our breakpoint analysis identifies multiple clouds in each region (Table \ref{tab:COMPLETE_distance}, 4th -- 6th columns).
On the other hand, \citet{2018ApJ...869...83Z} identifies distances for each velocity component of the $^{12}\mathrm{CO}~(1-0)$ emission line (Table \ref{tab:COMPLETE_distance}, 7th -- 13th columns).
Therefore, the clouds identified by this work and \citet{2018ApJ...869...83Z} do not necessarily have a one-to-one correspondence.
For simplicity, we show the closest pair of clouds by \citet{2018ApJ...869...83Z} in the same region for each cloud identified by the breakpoint analysis.
Similarly, for each cloud by \citet{2018ApJ...869...83Z}, the closest cloud in the same region by the breakpoint analysis is shown as a pair.
Therefore, note that the same cloud may be paired with two different clouds and shown twice in the figure.
}
\label{fig:COMPLETE_distance}
\end{figure}

\begin{table*}[t]
\begin{center}
\caption{Estimated Distances of Major Star-forming Regions across the Perseus Molecular Cloud}
\label{tab:COMPLETE_distance}
\begin{tabular}{lrr|rrr|rrrrrrr}
\hline
\hline
Cloud Name & R.A. & Decl. & \multicolumn{3}{c|}{This Work} & \multicolumn{7}{c}{\citet{2018ApJ...869...83Z}$^\mathrm{a}$}\\
& (deg) & (deg) & \multicolumn{3}{c|}{(pc)} & \multicolumn{7}{c}{(pc)} \\
\hline
L1451 & 51.0 & 30.5 & $146^{+39}_{-19}$ & ~ & $282^{+7~}_{-0~}$ & $66^{+3~}_{-2~}$ & $220^{+21}_{-16}$ & $223^{+11}_{-11}$ & $280^{+20}_{-22}$ & $292^{+6~}_{-4~}$ & & \\
L1448 & 51.2 & 30.9 & $154^{+61}_{-25}$ & $277^{+7~}_{-3~}$ & $290^{+7~}_{-3~}$ & $134^{+16}_{-40}$ & $286^{+9~}_{-8~}$ & $288^{+15}_{-17}$ & $289^{+10}_{-11}$ & $293^{+10}_{-11}$ & $294^{+11}_{-12}$ & \\
NGC 1333 & 52.2 & 31.2 & $176^{+41}_{-36}$ & $264^{+6~}_{-8~}$ & $294^{+34}_{-9~}$ & $110^{+24}_{-52}$ & $260^{+99}_{-44}$ & $287^{+20}_{-13}$ & $298^{+9~}_{-14}$ & $301^{+4~}_{-3~}$ & $302^{+5~}_{-6~}$ & $305^{+7}_{-7}$ \\
B1 & 53.4 & 31.1 & $142^{+41}_{-30}$ & $213^{+41}_{-30}$ & $300^{+6~}_{-1~}$ & $78^{+28}_{-13}$ & $266^{+28}_{-32}$ & $297^{+23}_{-33}$ & $300^{+4~}_{-4~}$ & $309^{+5~}_{-8~}$ & $313^{+6~}_{-5~}$ & \\
IC 348 & 55.8 & 31.8 & $131^{+15}_{-32}$ & $202^{+15}_{-32}$ & $309^{+8~}_{-2~}$ & $169^{+2}_{-3}$ & $294^{+5~}_{-4~}$ & $295^{+7~}_{-5~}$ & $295^{+5~}_{-4~}$ & $295^{+7~}_{-4~}$ & & \\
B5 & 56.9 & 32.9 & $156^{+10}_{-13}$ & ~ & $302^{+5~}_{-8~}$ & $79^{+16}_{-13}$ & $294^{+13}_{-11}$ & $302^{+8~}_{-9~}$ & $325^{+18}_{-17}$ & & & \\
\hline
\multicolumn{13}{p{0.85\textwidth}}{\bf Notes.}\\
\multicolumn{13}{p{0.85\textwidth}}{$^\mathrm{a}$Estimated distances to each velocity component of the CO molecular line emission for each cloud \citep{2018ApJ...869...83Z} is shown. The leftmost (closest) component of each cloud are foreground components.}\\
\end{tabular}
\end{center}
\end{table*}

We show the comparison of the results of our breakpoint analysis with the estimation by \citet{2018ApJ...869...83Z} in Figure \ref{fig:COMPLETE_distance} and Table \ref{tab:COMPLETE_distance}.
Our results are in reasonable agreement with those of \citet{2018ApJ...869...83Z} for each component of the Perseus cloud at a distance of about 300 pc.
The estimated distances show a gradual increase from west to east as pointed out by \citet[][also see \citealp{2019ApJ...877...88C}]{2018ApJ...869...83Z}.
Therefore, we can conclude that our \emph{G}-band breakpoint analysis can measure the distance to the cloud with sufficient accuracy.

On the other hand, the position of the foreground component is less consistent.
Therefore, we further check our estimated foreground distance's consistency, comparing that with 3D dust maps based on Gaia distance data by \citet{2019A&A...625A.135L} and \citet{2019A&A...631A..32L}.
We show a comparison between these results and ours for the Perseus molecular cloud direction in Figure \ref{fig:cloud_distance_profile}.

\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{cloud_distance_profile.png}
\caption{
Comparison of the distributions of the breakpoints distances in the Perseus molecular cloud direction ($157^\circ < l < 161^\circ,~-23^\circ < b < -16^\circ$) with existing 3D dust maps as a function of the distance.
Vertical axes are for the \emph{G}-band extinction ($A_G$) per each detected cloudlet (i.e., breakpoint) for panel (a), and $A_G$ per pc for panels (b) and (c).
}
\label{fig:cloud_distance_profile}
\end{figure}

Figure \ref{fig:cloud_distance_profile}a shows the results of our breakpoint analysis.
The spatial resolution of this analysis is set as $55\arcmin$.
The estimated cloud distances are concentrated at 150 pc and 300 pc (Figure \ref{fig:cloud_distance_profile}a).

Figures \ref{fig:cloud_distance_profile}b and c show the results by \citet{2019A&A...625A.135L} and \citet{2019A&A...631A..32L}, respectively.
Both two dust maps successfully detect the foreground cloud in addition to the Perseus cloud.
\citet[][Figure \ref{fig:cloud_distance_profile}b]{2019A&A...625A.135L} combine 2MASS photometric data with Gaia DR2.
Their map has limited LOS resolution of 50 pc and thus the estimated extinction for the Perseus cloud is smoothed out, resulting in relatively lower extinction per unit length value ($\sim 0.01$ mag/pc) compared to \citet[][Figure \ref{fig:cloud_distance_profile}c]{2019A&A...631A..32L}.
\citet{2019A&A...625A.135L} indicated plans to improve the LOS resolution as the Gaia data is updated.

\citet[][Figure \ref{fig:cloud_distance_profile}c, also see \citealp{2020A&A...639A.138L}]{2019A&A...631A..32L} used the Gaia DR2 $A_G$ values as photometric data to create a 3D dust map of the solar vicinity ($d \lesssim 300$ pc).
They achieved higher resolution in LOS comparing to the map by \citet{2019A&A...625A.135L} that used NIR photometric data (Figure \ref{fig:cloud_distance_profile}b).
NIR photometry is effective for tracing dense molecular clouds' interior but has limited sensitivity for tracing faint clouds with a high spatial resolution.
On the other hand, the disadvantage of tracing the shape of the dust cloud using only photometry of visible wavebands is also apparent; in their analysis shown in Figure \ref{fig:cloud_distance_profile}c, the 300 pc Perseus molecular cloud cannot be correctly detected due to saturation of the $A_G$ value and is split into two distances before and after the cloud.

Our results shown in Figure \ref{fig:cloud_distance_profile}a are based on the simple 2D sheet assumption and have limited spatial resolution of about $1^\circ$, but the use of the $A_G$ value achieves a good sensitivity to faint clouds as well as the correct determination of the distances of dense molecular clouds.

According to the discussions above, we judge that our method is a simple and effective one for measuring dust clouds.
Using \emph{G}-band photometry, we can obtain the distance to both faint and dense clouds stably and accurately.
Our results described in Section \ref{sec:breakpoints} (Figure \ref{fig:AngPbreak}) are thus consistent that the stellar polarization in the Perseus molecular cloud's direction originates in two isolated clouds at 150 pc and 300 pc, which are both detected by \citet{2019A&A...625A.135L} and \citet{2019A&A...631A..32L}.

\subsection{Spatial Distribution of Individual Clouds}
\label{sec: Taurus foreground}

Figure \ref{fig:TAP40} shows the spatial distribution of the dust cloud at 150 pc and 300 pc by \citet{2019A&A...631A..32L} as contours.
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{TAP50a.png}
\includegraphics[width=1.05\linewidth]{TAP50b.png}
\caption{
The spatial distribution of the dust clouds at 150 pc and 300 pc estimated by \citet{2019A&A...631A..32L}.
The blue contours are for the dust cloud integrated between 100 pc and 200 pc (i.e., the foreground cloud component), and the red contours are for the dust cloud integrated between 250 pc and 350 pc (i.e., the Perseus cloud itself).
The contour levels are [0.1, 0.3, 0.5, \ldots] mag in $A_G$.
{\bf Top panel}: the distribution in the Perseus region.
The display area of the figure is the same as that of Figures \ref{fig:Goodman} and \ref{fig:PerseusGaiadist}.
The background gray scale is the hydrogen column density ($N_\mathrm{H}$), and the black line segments are the \emph{B}-field orientation observed by Planck ($\theta_\mathrm{Planck}$; see Figure \ref{fig:Goodman}).
The spatial resolution of $\theta_\mathrm{Planck}$ is set as $10'$ in the top panel.
{\bf Bottom panel}: the distribution in the Taurus-Perseus region.
The yellow rectangle indicates the area shown in the top panel.
The yellow dashed circle shows the dust cavity outline identified in the 3D spatial distribution of dust clouds.
See Section \ref{sec:PerOB2} for more details on dust cavity identification.
The background gray scale is the Planck-observed 353 GHz continuum intensity.
The black line segments are $\theta_\mathrm{Planck}$, whose spatial resolution is set as $1^\circ$ in the bottom panel.
}
\label{fig:TAP40}
\end{figure}
To avoid the effect of LOS splitting of the dense molecular cloud, we show the spatial distribution of the integrated $A_G$ magnitude in the 100--200 pc range for the 150 pc component and in the 250--350 pc range for the 300 pc component, respectively.
In addition to the Perseus cloud (the top panel of Figure \ref{fig:TAP40}), we also show the Taurus-Perseus molecular cloud complex in the bottom panel of Figure \ref{fig:TAP40}.
As is evident in Figure \ref{fig:TAP40}, the 150 pc foreground cloud lying in front of the Perseus molecular cloud corresponds to the outer edge of the Taurus molecular cloud at a distance of $\sim 140$ pc \citep{2019ApJ...885...19Y, 2019ApJ...879..125Z, 2020A&A...638A..85R}, as was proposed by \citet{1987ApJS...63..645U} and \citet{1990Ap&SS.166..315C}.
The estimated extinction of the foreground cloud based on \citet{2019A&A...631A..32L} is $A_G^\mathrm{150pc} = 0.1$ -- 0.9 mag.

We overlay the \emph{B}-field orientation measured by Planck ($\theta_\mathrm{Planck}$) in Figure \ref{fig:TAP40} as black line segments.
$\theta_\mathrm{Planck}$ toward the Perseus cloud is generally aligned to the northwest-southeast direction.
The circular mean and the circular standard deviation are $\theta_\mathrm{Planck} = -56^\circ \pm 25^\circ$ if we estimate them  in the region where the 300 pc dust cloud component $A_G^\mathrm{300pc} > 0.6$ mag.
This $\theta_\mathrm{Planck}$ is consistent with the optical $\theta_\mathrm{star}$ of Group 1 ($-37.6^\circ \pm 35.2^\circ$).
Outside the Perseus molecular cloud but on the foreground component, we note a significant change in the \emph{B}-field orientation.
If we estimate $\theta_\mathrm{Planck}$ on the foreground cloud whose $A_G^\mathrm{150pc} > 0.2$ mag, $A_G^\mathrm{300pc} < 0.6$ mag, and R.A. $< 4^\mathrm{h}00^\mathrm{m}$, the position angle is $\theta_\mathrm{Planck} = -89^\circ \pm 19^\circ$, which is considerably different from $\theta_\mathrm{Planck}$ on the Perseus molecular cloud and is consistent with the $\theta_\mathrm{star}$ of Group 2 ($+66.8^\circ \pm 19.1^\circ$) within the range of errors.

As a result, we conclude that there are two \emph{B}-field components observed in the Perseus molecular cloud's direction: one is the Perseus molecular cloud's \emph{B}-field, and the other is that of a tenuous cloud at the outer edge of the Taurus molecular cloud with $A_G < 1$ mag.

\subsection{Polarization of the two clouds}
\label{sec:qu_separation}

\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{QU_separation.pdf}
\caption{
{\bf Left panel (a)}: the distribution of the observed relative Stokes $u$ and $q$ parameters, estimated from the observed $P_\mathrm{star}$ and $\theta_\mathrm{star}$ of Group 1 ($d > 300$ pc) and Group 2 ($150 < d < 300$ pc) stars.
The radial error bars show the observed error in $P_\mathrm{star}$, and the circumferential error bars show the observed error in $\theta_\mathrm{star}$.
The error bars may be hidden behind the marks.
{\bf Right panel (b)}: decomposition of the Group 1 $q$-$u$ vector into that of the Taurus cloud and the Perseus cloud.
The relationship between the average values of Taurus and Perseus clouds is shown.
}
\label{fig:QU_separation}
\end{figure}

The two-component \emph{B}-field of Perseus and Taurus clouds, described in the previous section, is traced by Group 1 ($d > 300$ pc) and Group 2 ($150 < d < 300$ pc) stellar polarization data.
Group 2 traces only the Taurus cloud, while Group 1 sees through both Perseus and Taurus clouds.
Here we estimate the Perseus molecular cloud's \emph{B}-field by removing the Taurus contribution from Group 1.

The observed polarization fraction is $\ll 10$ \% (Figure \ref{fig:AngPbreak} and Table \ref{tab:ANG_P_Group}).
In such a low polarization condition, the relative Stokes parameters $q$ (= Q/I) and $u$ (= U/I) in the polarized flux can be approximated as additive (e.g., \citealp{2019ApJ...872...56P} and the references therein).
In other words, we can separate the contribution to the polarization from individual clouds as follows.
\begin{eqnarray}
q_\mathrm{\,Group1} &=& q_\mathrm{\,Taurus} + q_\mathrm{\,Perseus}, \nonumber\\
q_\mathrm{\,Group2} &=& q_\mathrm{\,Taurus}, \nonumber
\end{eqnarray}
where $q_\mathrm{\,Group1}$ and $q_\mathrm{\,Group2}$ are the relative Stokes $q$ parameter observed for Group 1 and Group 2 stars, and $q_\mathrm{\,Taurus}$ and $q_\mathrm{\,Perseus}$ are the relative Stokes $q$ parameter originated in Taurus and Perseus clouds, respectively.
These formulas also hold for the relative Stokes $u$ parameter.
We estimate $q_\mathrm{\,star}$ and $u_\mathrm{\,star}$ values of Group 1 and Group 2 from the observed $\theta_\mathrm{star}$ and $P_\mathrm{star}$ as follows.
\begin{eqnarray}
q_\mathrm{\,star} &=& P_\mathrm{star} \times \cos\left(2 \cdot \theta_\mathrm{star}\right),\nonumber \\
u_\mathrm{\,star} &=& P_\mathrm{star} \times \sin\left(2 \cdot \theta_\mathrm{star}\right).\nonumber
\end{eqnarray}
The estimated $q_\mathrm{\,star}$ and $u_\mathrm{\,star}$ are shown in Figure \ref{fig:QU_separation}a.

The $q$ and $u$ values of Group 1 show a significant scatter, reflecting the local variation in $\theta_\mathrm{star}$ and $P_\mathrm{star}$ of Group 1.
On the other hand, the data of Group 2, which represents the $q$ and $u$ values of the Taurus cloud, show a small variation.
This is because $\theta_\mathrm{star}$ and $P_\mathrm{star}$ of Group 2 have a small spatial variation, and the value of $P_\mathrm{star}$ is also small (see Table \ref{tab:ANG_P_Group}).

Since $q$ and $u$ are additive, the observed values of $q$ and $u$ in Group 1 can be expressed as a vector sum of the contributions from Perseus and Taurus clouds on the $q$-$u$ plane.
This relationship is shown in Figure \ref{fig:QU_separation}b.
Here we show the relationship between the averaged $q$ and $u$ values of Taurus and Perseus.

The estimated $q$-$u$ vectors of Perseus and Taurus are nearly opposite to each other.
That is, the estimated orientations of the POS \emph{B}-field of these two clouds are nearly perpendicular to each other ($\theta_\mathrm{Perseus} \perp \theta_\mathrm{Taurus}$).
Furthermore, the Perseus contribution to the Group 1 $q$-$u$ vector is dominant compared to that of Taurus.
This may reflect the fact that the column density of the Perseus main cloud is dominant compared to that of the foreground Taurus cloud's outskirt (see Figure \ref{fig:cloud_distance_profile}).
As a result, the Taurus contribution to the Group 1 polarization does not significantly change the Perseus' q-u vector direction (Figure \ref{fig:QU_separation}b).
Thus, we conclude the followings.
\begin{eqnarray}
\theta_\mathrm{Perseus} &\simeq& \theta_\mathrm{Group1}, \nonumber\\ \theta_\mathrm{Taurus} &=& \theta_\mathrm{Group2}. \nonumber
\end{eqnarray}
On the other hand, we estimate that $P_\mathrm{Perseus}$ is slightly larger than $P_\mathrm{Group1}$ due to the depolarization by the foreground Taurus cloud.
We estimate $\theta_\mathrm{Perseus}$ and $P_\mathrm{Taurus}$ by subtracting the averaged Group 2 $q$-$u$ vector from the observed Group 1 $q$-$u$ vectors, as shown in Figure \ref{fig:QU_separation}b.
The results are shown in Table \ref{tab:p_efficiency}.

\begin{table}[t]
\begin{center}
\caption{Polarization properties of individual clouds}
\label{tab:p_efficiency}
\begin{tabular}{lcccc}
\hline
\hline
Cloud & $\theta_\mathrm{star}$ & $P_\mathrm{star}$ & $A_G$ & $P_\mathrm{star}/A_G$\\
& (deg) & (\%) & (mag) & (\%/mag) \\
\hline
Taurus & $+66.8 \pm 19.1$ & $0.8 \pm 0.6$ & $0.32 \pm 0.21$ & $1.5 \pm 0.3$\\
Perseus & $-30.0 \pm 25.2$ & $2.4 \pm 1.8$ & $1.61 \pm 0.58$ & $1.5 \pm 0.2$\\
\hline
\end{tabular}
\end{center}
\end{table}

Figure \ref{fig:pol_efficiency} shows $P_\mathrm{Perseus}$ and $P_\mathrm{Taurus}$ as a function of $A_G$ to investigate individual clouds' polarization efficiency ($P_\mathrm{star}/A_G$).
We assume that the Taurus cloud's $A_G$ values are equal to the Group 2 data.
For the Perseus cloud, we estimate $A_G$ by subtracting that of 150 pc cloud component estimated by \citet[][Figure \ref{fig:TAP40}]{2019A&A...631A..32L} at each position of the stars from the observed Group 1 $A_G$ values.
Note that the maximum observed $A_G$ value is 3 mag (also see Figure \ref{fig:offset_AG} and the discussion in Sections \ref{sec:PlanckCorrelation} and \ref{sec:small_scale}), indicating that the stellar extinction data are biased to regions in the cloud where the extinction is smaller ($A_G < 3$ mag).
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{pol_efficiency.pdf}
\caption{
$P_\mathrm{star}$ of Perseus and Taurus clouds as a function of $A_G$.
See text for the estimation of $P_\mathrm{star}$ and $A_G$ of individual clouds.
The dashed lines show the observed maximum polarization efficiency ($P/E(B-V) = 13 \%$; \citealp{2019A&A...624L...8P,2020A&A...641A..12P}) and its 50 \% slope.
We assume $A_G/A_V = 0.789$ and $A_V/E(B-V) = 3.16$ \citep{2019ApJ...877..116W} to convert $E(B-V)$ to $A_G$.
}
\label{fig:pol_efficiency}
\end{figure}

In Figure \ref{fig:pol_efficiency}, we also show the observed maximum polarization efficiency taken from the literature ($P/E(B-V) = 13 \%$;  \citealp{2019A&A...624L...8P,2020A&A...641A..12P}).
The polarization efficiency of Taurus and Perseus corresponds to less than $\sim 50\%$ of the maximum efficiency.
We do not find a significant difference between the polarization efficiencies of Taurus and Perseus clouds.

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

\subsection{\emph{B}-field Orientation and the Per OB2 Bubble}
\label{sec:PerOB2}

In Section \ref{sec: Taurus foreground}, we identified that the 150 pc foreground cloud lying in front of the Perseus molecular cloud is the outer edge of the Taurus molecular cloud.
Both Taurus and Perseus molecular clouds are thought to be on a common large H\small{I} shell's outer edge \citep{1974IAUS...60..115S,2006A&A...451..539S,2019A&A...623A..16S}.
In Figure \ref{fig:cloud_distance_profile}, we find a low-extinction region between these two clouds.
In the following, we estimate the spatial extent of this low-extinction region.

Figure \ref{fig:Perseus_chamber} shows the 3D distribution of the dust cloud in the region containing Taurus and Perseus.
\begin{figure*}[tp]
\centering
\includegraphics[width=\linewidth]{Perseus_chamber.png}
\caption{
The 3D distribution of dust clouds in the Taurus-Perseus region.
Coordinates are the heliocentric galactic cartesian coordinates;
the X-axis points to the galactic center, the Y-axis points to the direction of galactic rotation (galactic plane $l=90^\circ$), and the Z-axis points to the galactic north pole.
Points are the cloud positions determined by our breakpoint analysis, and the color scale is the distribution of dust clouds estimated by \citet[][available for the regions $\mathrm{X} \geq -300$ pc]{2019A&A...631A..32L}.
The distribution from parallel to the galactic plane at Z=0 pc to Z=-160 pc at every 20 pc is shown.
Yellow dashed ellipses show the cross-sectional profiles of the dust cavity (see Section \ref{sec:PerOB2}).
}
\label{fig:Perseus_chamber}
\end{figure*}
We estimate the 3D shape of the dust shell by using the cloud positions determined by our breakpoint analysis.
We assume the shell's profile as a simple ellipsoid and perform a least-square fit of the shell to the cloud positions weighted by their extinction (mag/cloud).
The estimated size and position of the shell are shown in Figure \ref{fig:Perseus_chamber_section}.
\begin{figure*}[tp]
\centering
\includegraphics[width=\linewidth]{Perseus_chamber_section.png}
\caption{
Cross-sectional views of the fitted ellipsoid to the dust shell are shown as dashed ellipses in panels (a): cross-sections in the direction of Y-axis, (b): same in X-axis, (c): same in Z-axis.
We show the cross sections in X, Y, and Z directions at the center of the fitted ellipsoid at $X = -220$ pc, $Y = 38$ pc, and $Z = -80$ pc.
The central position corresponds to $l = 170^\circ$, $b=-20^\circ$ in the galactic coordinates at the distance $d=240$ pc.
See section \ref{sec:PerOB2} for the description of the ellipsoid fit.
See Figure \ref{fig:Perseus_chamber} for the description of coordinates, color scales, and markers.
Dotted ellipses show the cross-sectional profiles of the dust cavity (see Section \ref{sec:PerOB2}).
The bottom right figures are radial profiles of the ISM, normalized to the radii of the fitted ellipsoid (the vertical dashed lines at radial distance = 1.0 correspond to the ellipsoid size in each direction).
Panel (d) is the distribution based on the breakpoint analysis, and panel (e) is that based on \citet{2019A&A...631A..32L}.
}
\label{fig:Perseus_chamber_section}
\end{figure*}
As shown in the figure, we can identify a low-density dust cavity surrounded by the dust shell.
The estimated diameters in the heliocentric galactic cartesian coordinates (X, Y, and Z directions; see Figure \ref{fig:Perseus_chamber} for the definition of X, Y, and Z) are
$D_X \simeq 156$ pc, $D_Y \simeq 110$ pc, and $D_Z \simeq 105$ pc, centering at $X \simeq -220$ pc, $Y \simeq 38$ pc, and $Z \simeq -80$ pc.
The ellipsoid's estimated shape is also shown in Figure \ref{fig:Perseus_chamber}, and the outline of the ellipsoid projected on the POS is shown in Figure \ref{fig:TAP40} (bottom panel).
The Taurus molecular cloud is located in front of this cavity, and the Perseus molecular cloud is located behind it.
Thus, we conclude that the two-component \emph{B}-fields observed in the Perseus cloud's direction at distances of 150 pc and 300 pc show the \emph{B}-field structure in front of and behind the dust cavity, respectively.

The H\small{I} shell and dust cavity are thought to be formed by the Per OB2 association.
The velocity field that is consistent with the assumed expanding motion of the shell was found for the Perseus molecular cloud in the H\small{I} \citep{1974IAUS...60..115S} and CO \citep{2006AJ....131.2921R,2006A&A...451..539S} line emissions (see, e.g., Figure 10 of \citealp{2006A&A...451..539S}).
\citet{2018A&A...614A.100T} found that the LOS \emph{B}-field directions are different from each other in the north and south of the east-west extending Perseus molecular cloud.
This LOS \emph{B}-field is consistent with a scenario where the environment has impacted and influenced the field lines to form a bow-shaped magnetic field morphology as explored in \citet{1997ApJS..111..245H} and \citet[][also see \citealp{2018PASJ...70S..53I}]{2019A&A...632A..68T}.
$\theta_\mathrm{star}$ of the Perseus cloud \emph{B}-field component is $-30.0^\circ \pm 25.2^\circ$ (Section \ref{sec:qu_separation} and Table \ref{tab:p_efficiency}).
In the galactic coordinates, $\theta_\mathrm{star}^\mathrm{GAL}$ (measured from Galactic North to Galactic East) of the Perseus cloud \emph{B}-field is $-68.5^\circ \pm 24.9^\circ$.
This $\theta_\mathrm{star}^\mathrm{GAL}$ is nearly parallel to the galactic plane (see Figure \ref{fig:AngPbreak}, top panel), i.e., it is consistent with the global \emph{B}-field component of the galaxy \citep[see, e.g.,][Figure 9]{2012ApJ...757...14J}. This $\theta_\mathrm{star}^\mathrm{GAL}$ is, therefore, in accordance with the formation scenario of the Perseus molecular cloud described above.
The impact of feedback on the \emph{B}-fields is explored on a forthcoming paper (Tahani et al., in prep).

On the other hand, the orientation of the Taurus \emph{B}-field is close to perpendicular to the galactic plane ($\theta_\mathrm{star}^\mathrm{GAL} = 28.7^\circ \pm 19.8^\circ$; Figure \ref{fig:AngPbreak}, top panel).
It is difficult to form a \emph{B}-field structure close to perpendicular to the galactic plane (i.e., perpendicular to the global \emph{B}-field of the galaxy) if the ISM is compressed simply along the LOS.
Thus, the observed $\theta_\mathrm{star}$ in front of and behind the dust cavity is not compatible with the simple expansion of H\small{I} bubble in a uniform \emph{B}-field.

As shown in the discussion in this section, thanks to the advent of Gaia data, we can now isolate the \emph{B}-fields associated with each of the multiple clouds superimposed along the LOS and map them to the 3D structure of the ISM \citep[see also][and the references therein]{2019ApJ...872...56P}.
This isolation of the \emph{B}-fields is important in studying the  \emph{B}-field properties of individual clouds, e.g., applying the DCF method \citep{1951PhRv...81..890D,1953ApJ...118..113C} to estimate the \emph{B}-field strength.

\subsection{Consistency between the Optical and the Planck-observed \emph{B}-field}
\label{sec:PlanckCorrelation}

In Figure \ref{fig:Goodman}, we showed the spatial distribution of $\theta_\mathrm{star}$ and $\theta_\mathrm{Planck}$.
We note that the Group 1 $\theta_\mathrm{star}$ ($d > 300$ pc; red points in Figure \ref{fig:Goodman}) has a good one-to-one correspondence with $\theta_\mathrm{Planck}$ in the same position, despite their spatial variations.
Figure \ref{fig:offset_distance} shows the offset angle between $\theta_\mathrm{star}$ and $\theta_\mathrm{Planck}$ as a function of stellar distance.
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{offset_distance}
\caption{
The angle difference between the optical polarimetry ($\theta_\mathrm{star}$) and the Planck-observed \emph{B}-field ($\theta_\mathrm{Planck}$) in the Perseus molecular cloud's direction.
The spatial resolution of the Planck observation is set as $10'$.
The offset angle is shown as a function of stellar distances observed by Gaia.
Note that the horizontal axis is a linear scale below 300 pc, but a logarithmic scale above 300 pc.
The dotted lines at $d = 150$ pc and 300 pc indicate the two breakpoints identified in $\theta_\mathrm{star}$ and $P$ (Section \ref{sec:breakpoints}), corresponding to the Taurus and Perseus clouds (Section \ref{sec: Taurus foreground}).
}
\label{fig:offset_distance}
\end{figure}
For stars closer than 300 pc, the offset angle is $\theta_\mathrm{star} - \theta_\mathrm{Planck} = -71^\circ \pm 45^\circ$ and shows large variation, while for stars farther than 300 pc, the offset angle is $\theta_\mathrm{star} - \theta_\mathrm{Planck} = 2^\circ \pm 28^\circ$.
$\theta_\mathrm{Planck}$ traces the \emph{B}-field approximately in proportion to the column density along the LOS.
Thus, $\theta_\mathrm{Planck}$ mainly traces the Perseus molecular cloud, with an additional contribution from the foreground cloud (see Figure \ref{fig:cloud_distance_profile}).
As discussed in Section \ref{sec:qu_separation}, this also applies to $\theta_\mathrm{star}$ of Group 1.
As a result, $\theta_\mathrm{Planck}$ agrees well with the Group 1 $\theta_\mathrm{star}$ as they trace the similar ISM on the LOS.
On the other hand, $\theta_\mathrm{star}$ of Group 2 and Group 3 has no contribution from the Perseus cloud, and the foreground contribution traced by them is much smaller than that of the Perseus cloud in the background.
Therefore, $\theta_\mathrm{star}$ of Group 2 and Group 3 does not correlate with $\theta_\mathrm{Planck}$.
Thus, we conclude that $\theta_\mathrm{star}$ and $\theta_\mathrm{Planck}$ observations of the Perseus molecular cloud's direction with stellar distances $d > 300$ pc are in good agreement with each other, despite the fact that their spatial resolutions are largely different from each other.
The spatial resolution of $\theta_\mathrm{Planck}$ is $10'$, which corresponds to $\sim 1$ pc at 300 pc from the sun.
On the other hand, the spatial resolution of $\theta_\mathrm{star}$ is the size of stellar photospheres, which is much smaller than the Planck's beam size.

The good one-to-one correspondence between $\theta_\mathrm{star}$ and $\theta_\mathrm{Planck}$ was first pointed out by \citet[][see also \citealp{2019ApJ...871L..15G}]{2015A&A...576A.106P}.
\citet{2016A&A...596A..93S} further analyzed the $\theta_\mathrm{star}$--$\theta_\mathrm{Planck}$ correlation and investigated possible small-scale structures of the \emph{B}-field inside the Planck beam.
They concluded that the \emph{B}-field inside the Planck beam is uniform or randomly fluctuated around the mean $\theta$, because $\theta_\mathrm{Planck}$ agrees with $\theta_\mathrm{star}$.
Our observed good correlation between $\theta_\mathrm{Planck}$ and $\theta_\mathrm{star}$ supports their conclusion.
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{offset_AG}
\caption{
The angle difference between the optical polarimetry ($\theta_\mathrm{star}$) and the Planck-observed \emph{B}-field ($\theta_\mathrm{Planck}$) in the Perseus molecular cloud's direction.
The spatial resolution of the Planck observation is set as $10'$.
The offset angle is shown as a function of stellar extinction $A_G$ observed by Gaia.
Corresponding $N_\mathrm{H}$ values are also shown at the top axis of the figure, assuming $N_\mathrm{H} = A_G \cdot 2.21\times 10^{21} / 0.789$ \citep{2009MNRAS.400.2050G,2019ApJ...877..116W}.
The colors of the points are the same as in Figure \ref{fig:AngPbreak}.
The dashed lines and the shaded area indicate the circular mean angle differences per $A_G = 0.5$ mag and their standard deviations.
}
\label{fig:offset_AG}
\end{figure}
We show the offset angle between $\theta_\mathrm{star}$ and $\theta_\mathrm{Planck}$ as a function of the stellar extinction in Figure \ref{fig:offset_AG}.
The data show good agreement (angle differences consistent with zero) with no clear dependence on column density.
The data shown here suggest that the small scale \emph{B}-field is smooth in the column density range up to $A_G \sim 3$ mag or hydrogen column density $N_\mathrm{H} \sim 10^{22}~(\mathrm{cm}^{-2})$.

In addition, Figure \ref{fig:offset_distance} suggests that the correlation between $\theta_\mathrm{Planck}$ and $\theta_\mathrm{star}$ doesn't change with the stellar distance beyond 300 pc.
As shown in Section \ref{sec:breakpoints} and Figure \ref{fig:AngPbreak}, $\theta_\mathrm{star}$ and $P_\mathrm{star}$ show no clear variation as a function of the stellar distance.
The good correlation between $\theta_\mathrm{star}$ and $\theta_\mathrm{Planck}$ thus indicate that the ISM between and behind the clouds does not significantly contribute to polarization in either emission ($\theta_\mathrm{Planck}$) or extinction ($\theta_\mathrm{star}$).
In other words, the polarization is produced only in dense clouds along the LOS.

\subsection{Small-scale \emph{B}-field in low- and high-column density regions}
\label{sec:small_scale}

\citet{2020ApJ...899...28D} performed sub-mm polarimetry observations of NGC 1333, a dense star-forming region in the Perseus molecular cloud.
They observed the region with POL-2 on JCMT, and revealed the \emph{B}-field orientation ($\theta_\mathrm{JCMT}$) with a high spatial resolution of $\simeq 0.02$ pc.
The observed region has ISM column density $N_\mathrm{H} \gg 10^{23}~\mathrm{(cm^{-2}}$; also see Figure \ref{fig:Goodman}).
In contrast to the good agreement between $\theta_\mathrm{star}$ and $\theta_\mathrm{Planck}$ described in the previous section, $\theta_\mathrm{JCMT}$ is not well correlated with $\theta_\mathrm{Planck}$.
Instead, $\theta_\mathrm{JCMT}$ shows significantly more complex structure on $< 0.5$ pc scales.
\citet{2020ApJ...899...28D} found that the small-scale \emph{B}-field appears to be twisted perpendicular to the gravitationally super-critical massive filaments, most probably due to the dense filaments' formation.

Hence, the small scale ($< 1$ pc) \emph{B}-field appears to be highly complicated for the regions with $N_\mathrm{H} \gg 10^{23}~(\mathrm{cm}^{-2})$ but smooth for the regions with $N_\mathrm{H} < 10^{22}~(\mathrm{cm}^{-2})$.
It is unknown about the small-scale \emph{B}-field structure for the regions with $N_\mathrm{H} = 10^{22}$ -- $10^{23}~(\mathrm{cm}^{-2})$.
This is because optical and NIR stellar polarimetry cannot trace $N_\mathrm{H} \gg 10^{22}~(\mathrm{cm}^{-2})$ ISM, while ground-based and airborne sub-mm telescopes currently lack sensitivity to trace $N_\mathrm{H} \ll 10^{24}~(\mathrm{cm}^{-2})$ ISM.
But as described below, observations of the LOS \emph{B}-field strength and Planck polarimetry with the lower spatial resolution hint at the increasing complexity of small-scale \emph{B}-field structure at $N_\mathrm{H} > 10^{22}~(\mathrm{cm}^{-2})$.
\citet{2012ARA&A..50...29C} showed that the LOS \emph{B}-field strength measured by Zeeman splitting measurements show a general increase with an increasing $N_\mathrm{H}$ when $N_\mathrm{H} > 10^{22}~(\mathrm{cm}^{-2})$.
In such high column density region, \citet{2015A&A...576A.104P} observed a sharp drop in $P_\mathrm{Planck}$ where $N_\mathrm{H} \geq 1.5 \times 10^{22}$ (cm$^{-2}$).
They also reported an increase of the polarization angle dispersion at the same column density range.
\citet{2016A&A...586A.138P} and \citet{2019A&A...629A..96S} found that the relative orientation of the Planck-observed \emph{B}-fields to the long axes of dense filaments changes systematically from being parallel to perpendicular at $N_\mathrm{H} \approx 5 \times 10^{21}$ (cm$^{-2}$), which is consistent with the result reported by \citet{2020ApJ...899...28D}.

Photometric observations also suggest that ISM clouds make a significant change in this column density range.
\citet{1998ApJ...502..296O} found recently formed protostars only in regions with $N_\mathrm{H_2} > 8 \times 10^{21}$ (cm$^{-2}$) in the Taurus molecular cloud.
\citet{2004ApJ...611L..45J} and \citet{2006ApJ...646.1009K} found no obvious substructures below an $N_\mathrm{H} \sim 10^{22}~(\mathrm{cm}^{-2})$ in their sub-millimeter continuum observations of the Ophiuchus cloud and the Perseus cloud, respectively.
Many authors have pointed out that there is a fiducial threshold of $N_\mathrm{H} \simeq 0.6$ -- $2 \times 10^{22}$ (cm$^{-2}$) for the star formation \citep[e.g.,][]{2010ApJ...724..687L,2012ApJ...745..190L,2010ApJ...723.1019H,2014ApJ...782..114E,2015A&A...584A..91K,2016MNRAS.459..342M,2019A&A...622A..52Z}, though some counterexample are reported by, e.g., \citet[][lower threshold value]{2020ApJ...904..172D} and \citet[][no threshold value]{2020ApJ...896...60P}.
Interestingly, the critical line mass of ISM filaments ($M_\mathrm{line,crit} = 2c_s^2/G \simeq 16~M_\odot~\mathrm{pc^{-1}}$ for $T_\mathrm{gas}=10$ K; \citealp{1963AcA....13...30S}; \citealp{1964ApJ...140.1056O}; \citealp{1997ApJ...480..681I}) is consistent with this threshold $N_\mathrm{H}$ value if we adopt the typical 0.1 pc width of the filaments \citep{2014prpl.conf...27A}.

These observations described above suggest that ISM clouds make a critical change above $N_\mathrm{H} \sim 10^{22}~(\mathrm{cm}^{-2})$.
Suppose the small scale \emph{B}-field becomes complex and the \emph{B}-field strength increases at $N_\mathrm{H} > 10^{22}~(\mathrm{cm}^{-2})$.
In that case, it could potentially show that the molecular clouds become gravitationally supercritical when $N_\mathrm{H} \gtrsim 10^{22}~(\mathrm{cm}^{-2})$, and thus the small scale structure in the molecular clouds are formed in this column density range.
The \emph{B}-field might be bent and distorted in the formation of gravitationally supercritical small scale structures or dense filaments.
Therefore, the \emph{B}-field structure may record these small-scale structures' formation history.
High spatial resolution (e.g., $< 15\arcsec$) interstellar \emph{B}-field observations in this column density range, if achieved, could provide important information on the formation of cloud structure in the early stages of star formation.

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

We studied the optical polarimetry data in the Perseus molecular cloud's direction together with the Gaia parallax distance of each star.
We found that observed values of both polarization angles and fractions show discrete jump at 150 pc and 300 pc and otherwise remain constant.
Thus, the polarization is originated in the dust clouds at 150 pc and 300 pc. The ISM between and behind the clouds does not make a significant contribution to the polarization.

The dust cloud at 300 pc corresponds to the Perseus molecular cloud.
We estimate the POS \emph{B}-field orientation angle of the Perseus molecular cloud as $-30.0^\circ \pm 25.2^\circ$.
The Perseus cloud has the highest column density in LOS, and thus Planck observations of the dust continuum in this direction mainly trace this cloud.
The optical polarization angles (stellar distances $d > 300$ pc) and the Planck-observed \emph{B}-field orientations are therefore well aligned with each other, although their spatial resolutions are largely different.

The dust cloud at a distance of 150 pc is faint with $A_G < 1$ mag, and the POS \emph{B}-field orientation angle is $+66.8^\circ \pm 19.1^\circ$.
We identified this cloud as the outer edge of the Taurus molecular cloud at the same distance.

The two dust clouds are at the front- and back-sides of a dust cavity, which corresponds to an H\small{I} bubble structure generally associated with the Per OB2 association.
We estimate the size of the dust cavity as 100 -- 160 pc.
The observed two components of the POS \emph{B}-field orientations show the \emph{B}-field orientations in front of and behind the cavity, which are nearly perpendicular to one another.

Thanks to the advent of Gaia data, we can now isolate the \emph{B}-fields associated with each of the multiple clouds superimposed along the LOS and map them to the 3D structure of the ISM.
This process can be applied to other regions, and we hope such efforts provide an important step toward understanding the 3D \emph{B}-field structure of the ISM.

\acknowledgments

This research has been supported by Grants-in-Aid for Scientific Research (25247016, 18H01250, 19H01938) from the Japan Society for the Promotion of Science.
This research was supported in part at the SOFIA Science Center, which is operated by the Universities Space Research Association under contract NNA17BF53C with the National Aeronautics and Space Administration.
M.M. is supported by JSPS KAKENHI grant No. 20K03276.
S.I.S acknowledges support from the Natural Science and Engineering Research Council of Canada (NSERC), RGPIN-2020-03981.
C.L.H.H. acknowledges the support of the NAOJ Fellowship and JSPS KAKENHI grants 18K13586 and 20K14527.
D.J. and R.P. are supported by the National Research Council  of  Canada  and  by  a  Natural  Sciences  and  Engineering  Research  Council  of  Canada  (NSERC)  Discovery Grant.
J.K.is supported by JSPS KAKENHI grant No. 19K14775.
M.T. is supported by JSPS KAKENHI grant Nos. 18H05442, 15H02063, and 22000005.
This work has made use of data from the European Space Agency (ESA) mission
{\it Gaia} (\url{https://www.cosmos.esa.int/gaia}), processed by the {\it Gaia}
Data Processing and Analysis Consortium (DPAC,
\url{https://www.cosmos.esa.int/web/gaia/dpac/consortium}). Funding for the DPAC
has been provided by national institutions, in particular the institutions
participating in the {\it Gaia} Multilateral Agreement.
This work has make use of observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada.
This research has also made use of the SIMBAD database and of NASA's Astrophysics Data System Bibliographic Services.

\facilities{Gaia, Planck}
\software{strucchange \citep{strucchange,breakpoints}, astropy \citep{2013A&A...558A..33A}}


\section{Introduction}

Quantum error correcting codes (or, in other terminology, quantum anticliques), introduced theoretically in the pioneering papers~\cite{shor, Ekert1996, Knill1997}, and implemented experimentally, e.g., in~\cite{Cory1998}, play an important role in quantum information theory~\cite{ManchiniWinterBook}. Analog error correction for continuous variables, such as position and momentum, was considered~\cite{Lloyd1998}, symmetry breaking in open quantum systems for photonic cat qubits~\cite{Lieu2020}, etc.  In particular, bosonic codes which use encoding of information in states of bosonis field are of high interest~\cite{Cochrane1999,Gottesman2001,Albert2018,TerhalQST2020,Albert2020}.

Error correction theory studies the possibility of encoding information in quantum states in a way to allow zero-error decoding in the presence of a given fixed acceptable set of errors. Mathematically each error is described by some completely positive map acting on the set of states of the quantum system. In the general setting~\cite{knill}, for any given set of errors it is possible to define a unique non-commutative operator graph \cite{Duan} such that the knowledge of this graph allows to define all error correcting codes for this set of errors. The correspondence between sets of errors and non-commutative operator graphs is not one-to-one, but in the case of separable Hilbert space each non-commutative operator graph describes codes for some set of errors~\cite{Duan2, shirokov, yashin}. As was established in the finite-dimensional \cite {Duan2, shirokov} and infinite-dimensional \cite {yashin} cases, the noncommutative graph describing errors in information transmission is always generated by some positive operator-valued measure.

Various models of error correction were analyzed using the approach based on the use of non-commutative graphs. It was applied for quantum error correction for the models of coupled finite-dimensional systems~\cite{AM1,AM2,AM3}, and for coupled infinite-dimensional system~\cite{AMP}. In \cite{AMP}, the non-commutative graph generated by the dynamics of a bipartite bosonic quantum system in an infinite dimensional Hilbert space was defined. The graph consists of orbits driven by the unitary group which is the solution of the Schr\"oedinger equation for a two interacting bosonic oscillators. In this framework possible error correcting codes are given by coherent states in the bosonic Fock space. In all these cases, it was possible to find quantum anticliques, which are projectors onto error correcting subspaces.

In this work, we extend the theory of operator graphs to the case when one system is finite dimensional while another is infinite dimensional. Currently, several quantum error-correction \cite{Kyaw1,Wu} and entanglement protection \cite{Fasihi} techniques was inroduced for systems of such structure. Explicitly, we consider the situation where the information is encoded in a joint state of a qubit (a two-level quantum system) coupled via the Jaynes-Cummings Hamiltonian~\cite{JayCum} with a bosonic oscillator or bosonic coherent field. The Jaynes-Cummings Hamiltonian is the key model for the various theoretical and experimental works in quantum optics and for studying the interactions between light and matter, see, e.g.~\cite{Rempe1987,Forn,Kockum,Casanova, Irish,Huang,Yoshihara,Langford,Cirac1993,Cirac1994,Haroche,Blais2004,Larson2007,Saffman2010}, including in strong~\cite{Cirac1993}, ultra-strong~\cite{Forn,Kockum,Huang} and deep strong coupling regimes~\cite{Forn,Kockum,Casanova,Huang,Yoshihara,Langford}.

We develop for this model the theory of non-commutative operator graphs and apply it to find the corresponding quantum anticlique. The construction is based upon Gazeau--Klauder coherent states~\cite{GeaKla}. With this setting, we explicitly find the error correcting subspace for any values of the parameters of the Jaynes-Cummings model.

The structure of the paper is the following. In section 2, we discuss the problem of finding the existence of an error correcting procedure for a quantum channel, and also describe quantum channels corresponding to operator graphs of the class which includes the non-commutative graphs later constructed in section 5. In section 3, the Jaynes-Cummings model is discussed. Section 4 describes the construction of Gazeau-Klauder coherent states. In section 5, using the Gazeau-Klauder coherent states for the Jaynes-Cummings model, we construct the non-commutative operator graphs that have quantum error correcting codes, and find quantum anticlique and error correction subspaces.

\section{Quantum channels and non-commutative graphs}

Consider encoding information in states of a quantum system with Hilbert space $\cal H$. The (convex) set  $\mathfrak {S}({\cal H})$ of quantum states is the set of positive unit trace operators in $\cal H$. Errors which can occur under information transmission can be described by a quantum channel $\Phi  :\mathfrak {S}({\cal H})\to \mathfrak {S}({\cal H})$ which is a completely positive trace preserving (CPTP) map. As any CPTP map, it possesses  the Kraus operator-sum representation (Kraus OSR) \cite {Hol}
\begin{equation}\label{count1}
\Phi (\rho )=\sum \limits _{k\in K}V_k\rho V_k^*,\quad \rho \in \mathfrak {S}(\mathcal{H}).
\end{equation}
The Kraus operators $\left\lbrace V_k,\ k\in K\right\rbrace $ are parametrized by some set $K$.
They should satisfy the property
\begin{equation}\label{count2}
\sum \limits _{k\in K}V_k^*V_k=\mathbb I.
\end{equation}
to preserve trace of density matrix. In infinite dimensional spaces $K$ is not necessary countable and the sum in (\ref {count1}) -- (\ref {count2}) can be replaced by an integral (see e.g. \cite{SI}). Nevertheless for any channel $\Phi $ there exists a countable set $K$ parametrizing Kraus operators such that (\ref {count1}) holds true. Note that the Kraus OSR of a quantum channel is non-unique. The same quantum channel also has a Kraus OSR with operators $\tilde V_j=\sum\limits_i U_{ji}V_i$, $j=1,\dots,m$ ($m\ge |K|$), where $U^+U=\mathbb I$ (unitary matrix).

The linear space, $\mathcal{V}$, plays an important role in the theory of optimal coding
\begin{equation}\label{gr}
\mathcal {V}=\overline {span} \left\lbrace V_k^*V_j,\ k,j\in K\right\rbrace .
\end{equation}
The linear space $\mathcal V$ does not depend on the choice of the operators $\{V_i\}$ used for Kraus OSR of a given quantum channel; despite the non-uniqueness of the Kraus OSR, it is unique for a given quantum channel. Notice that in the finite dimensional case there is not need to take the closure in~(\ref {gr}). For the infinite dimensional case see~\cite{AMP,yashin}.

The linear space  $\mathcal V$ has the properties of a non-commutative graph. Such objects were introduced in~\cite{ChoiEffros} as operator systems and recently redefined as non-commutative graphs in quantum information theory~\cite{Duan}.
A {\it non-commutative graph} is a linear subspace $\mathcal V$ of bounded operators in a Hilbert space ${\cal H}$ possessing the properties
\begin{itemize}
\item ${\bf V}\in \mathcal V$ implies that ${\bf V}^*\in \mathcal V$;
\item ${\bf I}\in {\mathcal V}$
\end{itemize}

The famous Knill--Laflamme condition\cite{Knill1997,knill} claims that a zero error transmission via some channel $\Phi $ is possible iff for some orthogonal projector $P$ for all $A\in \mathcal{V}$ holds $PAP=\alpha(A)P$, where $\alpha(A)\in \mathbb{C}$. Here $P$ is the projector on the subspace generated by error correction code~\cite{Knill1997}. The optimal code belongs to the subspace ${\cal H}_P=P{\cal H}$.  The dimension of the subspace ${\cal H}_P$ is the maximal amount of quantum information that could be transmitted via $\Phi$ with zero error. An orthogonal projection $P$ such that $\dim\left(  P {\cal H}\right) \ge 2$ is a {\it quantum anticlique} for a non-commutative graph $\mathcal V$ if it satisfies:
\begin{equation}\label {anti}
\dim P{\mathcal V} P=1.
\end{equation}

The most natural quantum channel is given by the projection measurement
\begin{equation}\label{channel}
\Phi _{\mathcal P}(\rho )=\sum \limits _{k\in K}P_k\rho P_k,\quad \rho \in \mathfrak {S}({\cal H}),
\end{equation}
where $\mathcal {P}=(P_k)$ is the orthogonal resolution of identity
$$
\sum \limits _{k\in K}P_k=\mathbb I.
$$
For the channel (\ref{channel}), the non-commutative graph (\ref{gr}) is
$$
\mathcal {V}=\overline {span} \left\lbrace P_k, k\in K \right\rbrace.
$$
In the finite dimensional case, it is enough to consider only discrete sets $K$, while for the infinite dimensional case $P_k=E(B_k)$ are generated by some projection valued measure on the real line, where $B_k\subset {\mathbb R}$ are some Borel sets possessing the property $\cup_{k\in K}B_k={\mathbb R}$. In this case, different choices of $B_k$'s produce different projectional measurements and a reachable set of admissible errors.

Suppose that some unitary group ${\mathcal U}=\{ U_t=e^{-itG},\ t\ge 0\}$ acts in the Hilbert space $H$. Then the possible set of errors can be extended to all  possible projection measurements
\begin{equation}\label{e*}
\Phi ^t(\rho )=\sum \limits _{k\in K}U_tP_kU_t^*\rho U_tP_kU_t^*,\ t\in {\mathbb R},
\end{equation}
For the goal  of constructing a quantum anticlique allowing to correct errors of the form (\ref {e*}) for any fixed $t$, it is natural to define the non-commutative graph corresponding to all these errors as follows
$$
{\mathcal V}=\overline {span}\{U_tP_kU_t^*,\ t\in {\mathbb R},\ k\in K\}.
$$
Our interpretation of such a graph is that the quantum system distorts the transmitted information by a set of time-dependent errors.

Based upon this interpretation, suppose that there is a set of orthogonal projections $\left\lbrace P_{\alpha },\ \alpha \in \mathfrak {A}\right\rbrace $ parameterized by some set $\mathfrak {A}$, and the operator space is generated by orbits of some unitary group $\mathcal U$ as follows
\begin{equation}\label{first}
{\mathcal V}=\overline {span}\{U_tP_{\alpha }U_t^*,\ t\in {\mathbb R},\ \alpha \in \mathfrak {A} \}.
\end{equation}
It is known \cite {Duan2, yashin} that (\ref {first}) is a non-commutative operator graph corresponding to some channel iff $\mathbb I\in \mathcal V$. We shall construct the explicit example of graph for the Jaynes-Cummings model. Moreover, we shall show that there exists an anticlique for this graph.

\section{Jaynes-Cummings model}

We consider a two-level quantum system (qubit) coupled to a coherent field. Hilbert space of the qubit is ${\cal H}_{\rm s}=\mathbb C^2$. Ground and excited basis states of the qubit are denoted as $\{\ket{g},\ket{e}\}$. Hilbert space of the coherent field is ${\cal H}_{\rm f}=L^2(\mathbb R)=\{f:\mathbb R\to\mathbb C\,|\, \int_\mathbb R|f(x)|^2 dx<\infty\}$. Fock states of the field are denoted as  $\{\ket{k}, k\in \mathbb{N}_0\}$. We use the set of natural numbers including zero $ \mathbb {N}_0 = {0} \cup \mathbb {N}$ to enumerate the states.
The qubit and the field are assumed to be coupled via the Jaynes-Cummings Hamiltonian acting in ${\cal H}={\cal H}_{\rm s}\otimes {\cal H}_{\rm f}$
\begin{equation}\label{HAM}
H=\omega_f a^{+}a^{-}+\frac{\omega_s}{2} \sigma_z + \frac{\kappa}{2}(\sigma^{-} a^{+}+\sigma^{+}a^{-}),
\end{equation}
Here $\omega_s, \omega_f \in \mathbb{R_+}$ are the frequencies of the qubit and the field, respectively, $\kappa\ge 0$ is the coupling constant, $\sigma_z$ is the Pauli matrix, $\sigma^{+},\sigma^{-}$ are the rising and lowering operators of the qubit and the $a^{+},a^{-}$ are the creation and annihilation operators of the field. The detuning parameter is $\Delta=\omega_f-\omega_s$. We use the normalization of the physical units such that $\hbar=1$. Denote the basis in ${\cal H}$ as $\ket{q}\otimes\ket{p}=\ket{q,p}$, where the first number $q\in \mathbb{N}_0$ denotes the coherent state and the second number $p\in \{e,g\}$ denotes the qubit state. The schematic picture of the Jaynes-Cumming model states interaction is provided in Fig.~\ref{figureJC}.

\begin{figure}[t]
\centering
\includegraphics[width=1\linewidth]{Fig1.png}
\caption{\label{figureJC} Jaynes-Cummings model of a qubit interacting with bosonic reservoir.}
\end{figure}

The Schr\"oedinger equation with the Jaynes-Cummings Hamiltonian has an exact solution. The Hamiltonian has the following eigenstates
\begin{eqnarray*}\label{dressedstates}
		&& \ket{0,g},\\
		&&\ket{n,+}=\cos \left(\frac{\theta_n}{2}\right)\ket{n-1,e}+\sin \left(\frac{\theta_n}{2}\right)\ket{n,g},\\
		&&\ket{n,-}=\sin \left(\frac{\theta_n}{2}\right)\ket{n-1,e}-\cos \left(\frac{\theta_n}{2}\right)\ket{n,g},
\end{eqnarray*}
where $\theta_n=\tan^{-1}(\kappa\sqrt{n}/\Delta)$ and $n\in \mathbb{N}$, for the non-resonant case $\Delta\neq 0$. For the resonant case $\Delta = 0$ the eigenstates are
\begin{eqnarray*}\label{dressedstates}
		&& \ket{0,g},\\
		&&\ket{n,+}=\ket{n-1,e}+\ket{n,g},\\
		&&\ket{n,-}=\ket{n,g}-\ket{n-1,e}.
\end{eqnarray*}
In both cases the corresponding eigenenergies are
\begin{eqnarray*}
		E_{0,g}&=& \frac{\omega_f+\Delta}{2} \\
		E_{n,\pm}&=&\omega_f\left(n-\frac{1}{2}\right)\pm \frac{1}{2}\sqrt{\Delta^2 + \kappa^2 n},\quad  n\in \mathbb{N}.
\end{eqnarray*}

Below we follow closely to the method provided in \cite{Hussin}, where a new class  of coherent states was constructed for the Jaynes-Cummings model with strictly increasing sequences of the eigenenergies $E_{n,+}$ and $E_{n,-}$. Our goal is to divide the space $\cal H$ into three direct summands, two of which are generated by eigenstates corresponding to strictly increasing sequences of eigenenergies and one is finite dimensional. The sequence
\begin{equation}\label {J}
		J_k=E_{k+1,+}, \quad k \in \mathbb{N}_0
\end{equation}
is known to be strictly increasing. On the other hand, the sequence $S_0=E_{0,g},S_k=E_{k,-},\ k\in\mathbb{N}$ may have degenerate levels $S_{k_1}=S_{k_2},\ k_1\neq k_2$. We want to keep only a strictly increasing tail of the sequence $S_k$. Let us show that there exists $M_0 \in \mathbb{N}$ such that for all  $l_2>l_1\ge M_0$ one gets $S_{l_2}> S_{l_1}$. It is equivalent to
\begin{widetext}
\begin{equation}
	S_{n+1}-S_n
	=\omega_f-\frac{1}{2}\left(\sqrt{\Delta^2+\kappa^2(n+1)}-\sqrt{\Delta^2+\kappa^2n}\right)>0,
	\quad \forall n\ge M_0.
\end{equation}
\end{widetext}
From this one gets that $M_0$ is the minimal integer solution of the inequality
\begin{equation}\label{K0}
		\left(\sqrt{\Delta^2+\kappa^2(M_0+1)}+\sqrt{\Delta^2+\kappa^2M_0}\right)^{-1}<\frac{2\omega_f}{\kappa^2}.
\end{equation}
Thus, the sequence
\begin{equation}\label{S}
		S_k=E_{k,-},\quad k\ge M_0,
\end{equation}
becomes strictly increasing.

Let us fix any number $K_0\in\mathbb N$, $K_0\ge M_0$. Then the sequence $ S_k=E_{k,-}$, $k\ge K_0$ will also be strictly increasing and we can  separate the pieces where we assured to have strictly increasing eigenenergies. This allows to represent the Hilbert space $\cal H$ as the direct sum ${\cal H}={\cal H}_1\oplus {\cal H}_2\oplus {\cal H}_3$, where
\begin{eqnarray}
		{\cal H}_1&=&span\{\ket{n,+},\ n\in \mathbb{N}\},\\
		{\cal H}_2&=&span\{\ket{n,-},\ n\ge K_0\},\\
		{\cal H}_3&=&span\{\ket{g,0}\}\cup\{\ket{n,-},\ 1\le n < K_0\},
\end{eqnarray}
The subspace ${\cal H}_3$ later will be shown to be the error correcting subspace for this system. The error dynamics will be shown to interchange states in ${\cal H}_1$ and ${\cal H}_2$, while keeping states in $\mathcal{H}_3$ unchanged.

In Section 5, following the ideas of~\cite{Hussin} we will define the Gazeau-Klauder coherent states in ${\cal H}_1$ and ${\cal H}_2$.

\section{Gazeau-Klauder coherent states}

Here we introduce the construction of Gazeau-Klauder coherent states \cite{GeaKla}. Let us consider an infinite dimensional Hilbert space ${\cal H}$ with the basis $\ket{k},k\in \mathbb{N}_0$ and a self-ajoint operator $G$ which is diagonal in this basis. In \cite{GeaKla} Gazeau and Klauder defined the generalized coherent states corresponding to the operator $G$ as a two-parameter system of vectors $\{\ket{x,y}, x\in \mathbb{R_+},\ y\in \mathbb{R} \}\subset \mathcal{H}$ with the following properties
\begin{enumerate}\label{gcs}
		\item \textit{Continuity:}$\ (x,y)\rightarrow(x_0,y_0)\Rightarrow\ket{x,y}\rightarrow \ket{x_0,y_0}$.
		\item \textit{Resolution of identity:} $\int \ket{x,y}\bra{x,y} d\nu(x,y)=\mathbb I_{\cal H}$.
		\item \textit{Temporal stability:} $e^{-itG}\ket{x,y}=\ket{x,y+\omega t}$.
		\item \textit{Action identity:} $\bra{x,y} G\ket{x,y}=\omega x$.
\end{enumerate}
for some real constant $\omega$ and some measure $\nu$.

Consider the set of eigenvalues $h_k=\bra{k}G\ket{k}$ for the operator $G$. In the case $h_0=0$ and strictly increasing $h_k$, Gazeau and Klauder gave the explicit construction for  the system of coherent states. If $h_0> 0$ and the sequence $h_k$ is strictly increasing, their construction describes the set of vectors that satisfy the first two properties and the following version of the time stability condition:
$$
e^{-itG}\ket{x,y}\bra{x,y}e^{itG}=\ket{x,y+\omega t}\bra{x,y+\omega t}.
$$

The main property of generalized coherent states is that they form the resolution of identity. It was shown that the measure $\nu $ has the form
$$
d\nu (x,y)=\tau (x)dxdy,
$$
where $\tau (x)$ is some probability distribution density on the half-axis.

Below we give an explicit description of this construction. Consider a sequences of weights
$$
c_k>0,\quad k \in \mathbb{N}_0,
$$
with the convergence condition
\begin{equation}\label{conv}
\limsup_{k\rightarrow \infty} \sqrt [k]{c_k}=R >0.
\end{equation}
Suppose that these weights are the moments of probability distributions with  density $\rho(x)> 0$ on the interval $[0,R)$,
$$
c_k=\int\limits_0^R \rho(x) x^k dx < +\infty,\quad  k \in \mathbb{N}_0.
$$
We also need the normalization factor and the density defined by the formulae
\begin{eqnarray}\label{ser}
N^2(x)&=&\sum_{k=0}^{\infty}\frac{x^k}{c_k},\quad  0\le x< R,\\
\tau(x)&=&N^2(x)\rho (x).\nonumber
\end{eqnarray}
The radius of convergence in (\ref {ser}) is equal to $R$ by the property (\ref{conv}).  Now the Gazeau--Klauder coherent states are defined as follows
\begin{equation}\label{GKS}
\ket{x,y}=\frac{1}{N(x)}\sum_{k=0}^{\infty}\frac{x^{k/2}e^{-ih_{k}y}}{\sqrt{c_k}}\ket{k}
\end{equation}
We suppose the constant $\omega$ is equal to one.
Following the definition in \cite{GeaKla}, for $f:\mathbb{R}\rightarrow B(\mathcal{H})$, where $B(\mathcal{H})$ is the set of bounded linear operators on $\mathcal{H}$, we introduce its integration as the weak-limit of averages of weak integrals,
$$
I(f)=\int\limits_{-\infty}^{+\infty}f(y)d\mu(y)=\lim_{R \rightarrow +\infty} \frac{1}{2R}\int\limits_{-R}^{R}f(y)dy.
$$
Note that if such integral converges for some $f$ with the image lying in the weakly-closed subspace $Im(f)\subset\mathcal{W}\subset B(H)$, then the integral also lies in this subspace, $I(f)\in\mathcal{W}$. The resolution of identity property for coherent states results in
\begin{equation}\label{resol}
\int\limits_{0}^{R}\int\limits_{-\infty}^{+\infty}\ket{x,y}\bra{x,y}\tau(x)dxd\mu(y)=\mathbb I_{\cal H}.
\end{equation}

\section{Graphs generated by Gezeau-Klauder coherent states}

Now we are able to define systems of Gazeau--Klauder coherent states \cite{GeaKla} in the Hilbert spaces ${\cal H}_1$ and ${\cal H}_2$. Take two sequences of weights
$$
c_k^{(j)}>0,\quad d \in \mathbb{N}_0,\   j=1,2,
$$
with the same convergence condition
$$
\limsup_{k\rightarrow \infty} \sqrt[k]{c_k^{(j)}}=R >0.
$$
Suppose that the weights have the corresponding probability densities $\rho_1(x),\rho_2(x)> 0$ on the interval $[0,R]$ such that
$$
c_k^{(j)}=\int\limits_0^R \rho_j(x) x^k dx < +\infty,\quad k \in \mathbb{N}_0,\   j=1,2.
$$
Then, the normalization factors and the densities for measures defining resolutions of identity are given by the formulae
\begin{eqnarray*}
N_{j}^2(x)&=&\sum_{k=0}^{\infty}\frac{x^k}{c_k^{(j)}}, \quad 0\le x <R,\\
		\tau_{j}(x)&=&N_{j}^2(x)\rho_j(x),\quad j=1,2.
\end{eqnarray*}
Consider the Gazeau--Klauder coherent states
\begin{eqnarray*}
		\ket{J,x,y}&=&\frac{1}{N_1(x)}\sum_{k=0}^{\infty}\frac{x^{k/2}e^{-iJ_ky}}{\sqrt{c_k^{(1)}}}\ket{k+1,+},\\
		\ket{S,x,y}&=&\frac{1}{N_2(x)}\sum_{k=0}^{\infty}\frac{x^{k/2}e^{-iS_{k+K_0}y}}{\sqrt{c_k^{(2)}}}\ket{k+K_0,-},
\end{eqnarray*}
where the strictly increasing sequences of eigenenergies are given by (\ref {J}) and (\ref {S}) respectively.

Since the Gazeau--Klauder coherent states section \ref {gcs} form the resolution of identity, we get for the projections $P_{{\cal H}_1},P_{{\cal H}_2}$ on the subspaces ${\cal H}_1, {\cal H}_2$ that
\begin{equation}\label{resol1}
		\int\limits_{0}^{R}\int\limits_{-\infty}^{+\infty}\ket{J,x,y}\bra{J,x,y}\tau_1(x)dx d\mu(y)=P_{{\cal H}_1},
\end{equation}
\begin{equation}\label{resol2}
		\int\limits_{0}^{R}\int\limits_{-\infty}^{+\infty}\ket{S,x,y}\bra{S,x,y}\tau_2(x)dxd\mu(y)=P_{{\cal H}_2}.
\end{equation}

Consider the unitary group ${\mathcal U}=\{ U_t=e^{-it{H}},\ t\in {\mathbb R}\}$, where the Hamiltonian $H$ is determined by~(\ref {HAM}). Systems $\ket{J,x,y}, \ket{S,x,y}$ satisfy the temporal stability property~(\ref {gcs}) with respect to $\mathcal U$,
\begin{equation}\label{ts1}
		U_t\ket{J,x,y}\bra{J,x,y}U_t^{*}=\ket{J,x,y+t}\bra{J,x,y+t}
\end{equation}
\begin{equation}\label{ts2}
		U_t\ket{S,x,y}\bra{S,x,y}U_t^{*}=\ket{S,x,y+t}\bra{S,x,y+t}
\end{equation}

Consider the two families of orthogonal projections
$$
P_x^1=\ket{J,x,0}\bra{J,x,0},\quad  P_x^2=\ket{S,x,0}\bra{S,x,0},
$$
for $x\in [0,R]$. The projections $P_x^1$, $P_x^2$ and $P_x^3\equiv P_{{\cal H}_3}$ are pairwise orthogonal for any fixed value of $x\in [0,R]$.

{\bf Theorem 1.} The subspace
$$
	\mathcal{V}=\overline{span}\lbrace U_tP_x^jU_t^{*},\ t\in \mathbb{R}, \ x \in [0,R] ,\ j\in \{1,2,3\}\rbrace
$$
is a non-commutative operator graph with the anticlique $P_{{\cal H}_3}.$

\noindent {\bf Proof.} Consider the operator
$$
Q_x=P_x^1+\frac {\tau _2(x)}{\tau _1(x)}P_x^2+\frac {1}{\tau _1(x)}P_{{\cal H}_3}
$$
It follows from (\ref {ts1}) and (\ref {ts2}) that
\begin{eqnarray*}
	U_tQ_xU_t^{*}\\
	=\ket{J,x,t}\bra{J,x,t}+\frac{\tau_2(x)}{\tau_1(x)}\ket{S,x,t}\bra{S,x,t}+ \frac{1}{\tau_1(x)}P_{{\cal H}_3}.
\end{eqnarray*}
Then, (\ref{resol1}) and (\ref {resol2}) result in
\begin{widetext}
\begin{equation}
				\int\limits_{0}^{R}\int\limits_{-\infty}^{+\infty}\tau_1(x)\left(\ket{J,x,t}\bra{J,x,t}
			+ \frac{1}{\tau_1(x)}\mathbb I_{{\cal H}_2}+ \frac{\tau_2(x)}{\tau_1(x)}\ket{S,x,t}\bra{S,x,t}\right)
			dxd\mu(t)=\mathbb I_{\cal H}\in {\mathcal V}.
\end{equation}
\end{widetext}
Since $K_0$ is given by the rule (\ref{nontriv}) the dimension of ${\cal H}_3$ is at least 2. From the equalities
\begin{eqnarray*}
		P_{{\cal H}_3}\ket{J,x,t}\bra{J,x,t}P_{{\cal H}_3}&=&0\\
		P_{{\cal H}_3}\ket{S,x,t}\bra{J,x,t}P_{{\cal H}_3}&=&0
\end{eqnarray*}
we obtain that $P_{{\cal H}_3}$ is an anticlique. $\Box$

\section{The error correcting subspace}
As Theorem~1 states, the subspace ${\cal H}_3$ is the error correcting subspace. In some cases the number is $M_0=1$, so in this case for $K_0=M_0$ the error correcting subspace would be empty (since its dimension is $K_0-1$). However, in our construction one can take any natural $K_0\ge M_0$. To satisfy this condition for our coding procedure, that is the dimension of error correcting subspace is greater or equal to two, we should take any

\begin{equation}\label{nontriv}
		K_0\ge K_0^*=\max\{3,M_0\}.
\end{equation}

To analyze the minimal dimension of the error correcting subspace for various parameters of the Jaynes-Cummings Hamiltonian, consider the coupling rates $\gamma_f=\kappa/ \omega_f$ and $\gamma_s=\kappa/ \omega_s$. In terms of these quantities the inequality~(\ref{K0}) takes the following form
\begin{widetext}
\begin{equation} \label{mindim}
	 \sqrt{(\gamma_f^{-1}-\gamma_s^{-1})^2+M_0+1}+\sqrt{(\gamma_f^{-1}-\gamma_s^{-1})^2+M_0} >\frac{\gamma_f}{2}.
\end{equation}
\end{widetext}
Now it is evident that the key quantity defining the dimension for the resonant case (when  $\Delta=\omega_f-\omega_s=0$ and hence $\gamma_f=\gamma_s$) is the coupling rate $\gamma_f$. For having $M_0$ equal or larger than $4$ or greater, we need the $\gamma_f$ to be at least $2(2+\sqrt{3})$. For non-resonant case for fixed $\gamma_f$ decreasing the value $\gamma_s$ will decrease $M_0$. Figure \ref{figureBeh} shows the behavior of the minimal possible dimension $D_{\rm min}=K_0^*-1$ of the error correcting subspace ${\cal H}_3$ vs coupling rates $\gamma_s$ and $\gamma_f$. The figure clearly shows that the behavior is non-symmetric with respect to $\gamma_s$ and $\gamma_f$, as is also evident from inequality~(\ref{mindim}). The resonant case is the extremal case in the inequality~(\ref{mindim}), what means if $M_0$ satisfies it for some $\gamma_f$, then for the same parameter $\gamma_f$ the number $M_0$ also solves the inequality in the resonant case. Figure~\ref{figureRes} shows the dependence of $D_{\rm min}$ on the coupling rate $\gamma_f$ for equal frequencies.

\begin{figure}[h]
		\includegraphics[width=1\linewidth]{Fig2.png}
		\caption{\label{figureBeh} Behavior of the minimal dimension of the error correcting subspace vs coupling rates of the Jaynes-Cummings Hamiltonian.}
\end{figure}

\begin{figure}[h]
		\includegraphics[width=1\linewidth]{Fig3.png}
		\caption{\label{figureRes} Behavior of the minimal dimension of the error correcting subspace for the resonant case $\Delta=0$ vs coupling rate of the field.}
\end{figure}

The Jaynes-Cummings Hamiltonian is used in various theoretical and experimental analysis in quantum optics, cavity QED, e.g.~\cite{Rempe1987,Forn,Kockum,Casanova, Irish,Huang,Yoshihara, Langford,Cirac1993,Cirac1994,Haroche,Blais2004,Larson2007,Saffman2010}, including in strong~\cite{Cirac1993}, ultra-strong~\cite{Forn,Kockum,Huang} and deep strong coupling regimes~\cite{Forn,Kockum,Casanova,Huang,Yoshihara,Langford}. We take for example of the weak coupling regime the parameters used in experiments performed in the group of Serge Haroche~\cite{Haroche}. In the setup of that experiments the cavity is designed to have the frequency equal to the frequency of the atom, i.e. $\Delta = 0$.   The approximate parameters for the experiment are $\kappa=2\pi \cdot 47$ KHz and $\omega_f=\omega_s=2\pi\cdot 51.1$ GHz. For this case the inequality~(\ref{K0}) becomes
\begin{equation} \label{mindim1}
	\sqrt{M_0+1}+\sqrt{M_0}>\frac{\kappa}{2\omega_f}=\frac{47}{2\cdot51.1\cdot10^6} \approx 0.46 \cdot 10^{-6}
\end{equation}
The minimal natural solution of this inequality is $M_0=1$, so $K_0=3$ and the minimal dimension is $D_{\rm min}=K_0-1=2$. This minimal two-dimensional error correcting subspace is spanned by the two vectors $\ket{g,0}$ and $\ket{1,-}=\ket{1,g}-\ket{0,e}$. We remark that this is the error correcting subspace of minimal dimension. One could choose arbitrary large $K_0$ and the corresponding non-commutative operator graph will have the error correcting subspace $H_3$ of dimension $K_0-1$.

Jaynes-Cummings model is derived from Rabi model via the rotating wave approximation (RWA). This approximation is typically valid for $\gamma_f<0.1$ and $\omega_f\approx \omega_s$. In this case the minimal dimension is $D_{\rm min}=2$.  As it can be seen from Fig.~\ref{figureBeh} and  Fig.~\ref{figureRes}, to have $D_{\rm min}>2$ one has to consider values of $\gamma_f$ in the range of deep strong coupling regime~\cite{Forn,Kockum,Casanova}. This regime, as well as less intense ultrastrong regime, is of interest now. In these regimes the Rabi model is non-integrable and investigating these regimes motivates describing eigen-energies approximations for this model~\cite{Casanova,Irish}. One can show \cite{Huang} that introducing a special type of frequency modulations applied to the field and the qubit will give the dynamics governed by Jaynes-Cummings Hamiltonian with $\gamma_f$ in the range of deep strong coupling regime. In circuit-QED simulations rates of $\gamma_f$ for Rabi model up to 2.1 are achieved~\cite{Yoshihara, Langford}. Thus value is lower, while not that much, than the minimal $\gamma_f\approx 7.5$ that is necessary to see the effect in which minimal dimension of the error correcting code in the proposed scheme will be 3 or greater. Our analysis  allows to construct non-trivial quantum error correcting codes for all possible values $\gamma_f,\gamma_s$.

We remark that our scheme could be applied to any system that possess the same decomposition into the direct sum, where eigenenergies in the two direct summands form strictly increasing sequences. Potentially, this property could be exploited for more complex Hamiltonian beyound the Jaynes-Cummings model, as for example for Jaynes-Cummings-Hubbard Hamiltonian, describing interaction of several qubit-cavity systems, or for perhaps directly for Rabi Hamiltonian.

\section{Conclusion}
In this work, the theory of non-commutative operator graphs has been developed for error correction in the case of a finite-dimensional quantum system coupled to an infinite-dimensional quantum system. We have constructed the non-commutative operator graph generated by orbits of the unitary group driven by Hamiltonian (\ref {HAM}) of the Jaynes-Cummings model. We have shown that for a positive integer $K_0$ that satisfies (\ref{nontriv}), using for encoding the eigenstates $\ket {g,0}$ with $\ket {n ,-}$ for $0\le n<K_0$ (\ref {dressedstates}) allows to transmit information with zero error via quantum channels with operator graphs belonging to the constructed graph. Thus the error correcting subspace is explicitly computed for all values of the parameters of the Jaynes-Cummings model. Our scheme could be applied to any system that possess the same decomposition of eigenenergies into the direct sum as for JC Hamiltonian, where eigenenergies in the two direct summands form strictly increasing sequences.

\begin{acknowledgments}
This work was funded by the Ministry of Science and Higher Education of the Russian Federation (grant number 075-15-2020-788) and performed at the Steklov Mathematical Institute of the Russian Academy of Sciences.
\end{acknowledgments}