Exactness of the mean-field dynamics in optical cavity systems

Quantum dynamics in many body systems is a very important topic for many branches of physics. In particular, recent studies in the field of statistical physics have focused on dynamical problems such as thermalization [1]–[5], phase transitions induced by the parameter quench [6]–[9], and nonequilibrium phase transitions [10]–[12]. Theoretically, it is difficult to precisely describe the time evolution of a many body system. Although recent works have succeeded in treating the exact time evolution and observing relaxation processes in some integrable systems [13]–[15], we should be content with some approximate treatments in general.

One of the important approximate approaches is the mean-field (MF) approximation for not only equilibrium statics but also nonequilibrium dynamics. In this approximation, the N-body state vector or the N-body density matrix is approximated by the product state. Interestingly, the MF approach predicts some remarkable results such as the absence of thermalization and purely dynamical phase transitions in some isolated quantum systems [9, 10]. In open quantum systems, nonequilibrium phase transitions have also been studied by the MF approach [16]–[18]; see also [19] for a recent result.

Some of these predictions might be the peculiarity of the MF theory and not universal in general. However, importantly, there are several models in which the MF treatment becomes exact in some ideal limit [20]–[23]. Thus the MF theory is completely reliable as long as such an ideal limit is considered to be realized. For example, in a quantum spin system with a global coupling (the infinite-range interaction), it is known that the quantum dynamics is exactly described by the MF dynamical equation (or Hartree equation) [20]. Another example is an N-body bosonic system interacting via a two-body potential with scattering length a, in which the time-dependent Gross–Pitaevskii equation, which is regarded as the MF theory of the dynamics of the Bose–Einstein condensate, becomes exact in the limit of N → ∞ with a fixed value of Na [21]. Because of the existence of such ideal limits, the MF approximation takes a special position among other approximate approaches.

Therefore, we should not regard the MF approach as a mere crude approximation. Rather, it is important to clarify the conditions under which the MF approach is justified and extend the possibility of experimental realization of the remarkable predictions of the MF approach.

This paper is mainly devoted to investigating the quantum dynamics of a system consisting of N two-level atoms interacting with a single quantized boson mode in an optical cavity. This system is described by the Dicke model [24]. We consider its extension, that is, we allow each element of the system, i.e. the single cavity mode of the photons and N two-level atoms, to be attached to its own environment. Each environment is a Hamilton system, and may or may not be a large system. If one environment has a continuous energy spectrum, which means that the environment is an infinitely large system, it plays the role of a dissipative thermal bath.

The MF approach of this model predicts interesting dynamical or nonequilibrium phase transitions. In the isolated Dicke model (there is no environment), it was argued that dynamical phase transitions occur by the MF approach [9]. In open systems (the Dicke model interacting with dissipative environments), under the time-periodic driving field, the MF dynamics combined with the method of the Born–Markov quantum master equation [25] predicts some nonequilibrium phase transitions including optical bistabilities [17, 18], which have been observed in experiments [26]. It has been believed that the MF approach is justified in the limit of N → ∞ because in this model N two-level atoms interact with a common cavity mode of photons, and it is expected that each atom feels only the MF produced by the other atoms via radiation and absorption of cavity photons. In this situation, the correlation between two atoms or between an atom and the cavity photons will be not so important. However, there has not been rigorous proof of this expectation yet.

In this paper, we rigorously prove that the MF approach becomes exact in the limit of an infinite number of two-level atoms in the sense that the expectation value of an observable belonging to the restricted set ℬ, which will be specified in section 4.2, at an arbitrary time t is exactly equal to that calculated by the MF theory. The system may be in contact with thermal reservoirs and may be subject to a time-dependent external field, and the result does not depend on whether the dynamics of the system of interest is Markovian or not. Thus the result is very general. The use of the MF theory in the above mentioned previous works on this model is justified. The restriction of this work, which should be removed in future works, is that the Hamiltonian must be linear with respect to cavity photons and the Bose field attaching to them (the '0-subsystem' in section 3).

The strategy of the proof is based on comparison of the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) hierarchy generated by the exact equations with the one generated by the MF equations. We will see that they become indistinguishable in the thermodynamic limit. This strategy itself is applicable to simpler models, i.e. infinite-range spin models, so that we will first apply this strategy to those simple models, and then proceed to the extended Dicke model, in which the proof is much more complicated.

The organization of this paper is as follows. In section 2, we introduce the proof of the fact that the MF theory is justified in infinite-range spin models. This result is already known [20], but we review it because it is the simplest case where the MF approach is justified rigorously. In section 3, the model considered in this work is explained. In section 4, some preliminaries are covered. We introduce the two notions that are important later, i.e. the Q-representation of the density matrix and the restriction of the class of physical observables. In section 5, the main result of this paper is stated, and its proof is summarized. The detailed evaluation of some quantities necessary for the proof is given in appendices A and B. In section 6, we summarize the result of this work and discuss some future problems.

Before studying the optical cavity system, let us analyze the simplest case where the MF dynamics becomes exact in the thermodynamic limit. We consider an N spin system with infinitely long-range interactions. The Hilbert space for the ith spin is denoted by _i. Let h_i = h be a bounded self-adjoint operator acting on _i, and let V_ij = V (i ≠ j) also be a bounded self-adjoint operator acting on _i⊗_j. We put V_ii = 0. The Hamiltonian of the spin system with infinite-range two-body interactions is generally written as

$$\begin{eqnarray} &&H=\sum _{i=1}^{N}{h}_{i}+\frac{1}{N}\sum _{i,j=1}^{N}{V}_{i j}. \end{eqnarray} \tag{ 1 }$$

We introduce the p-norm of an operator A acting on ${\mathop{\bigotimes }\nolimits }_{i=1}^{k}{\mathcal{H}}_{i}$ as

$$\begin{eqnarray} &&\Vert A\Vert _{p}^{(k)}:= \left ({\mathrm{Tr}}_{1,2,\ldots ,k}\vert A\vert ^{p}\right )^{1/p}. \end{eqnarray} \tag{ 2 }$$

In particular, the operator norm is defined as

$$\begin{eqnarray} &&\Vert A\Vert _{\infty }:= \mathop{\sup }\limits _{\psi \in \bigotimes _{i=1}^{k}{\mathcal{H}}_{i}}\frac{\langle \psi \Vert A\Vert \psi \rangle }{\langle \psi \vert \psi \rangle }. \end{eqnarray} \tag{ 3 }$$

The operator A is said to be bounded if ||A||_∞ <+ ∞.

The N-spin density matrix at time t is denoted by ρ_N,t, and its k-marginal reduced density matrix by ${\gamma }_{N,t}^{(k)}$, which is defined as

$$\begin{eqnarray} &&{\gamma }_{N,t}^{(k)}:= {\mathrm{Tr}}_{k+1,k+2,\ldots ,N}{\rho }_{N,t}. \end{eqnarray} \tag{ 4 }$$

The initial condition is assumed to satisfy

$$\begin{eqnarray} &&\mathop{\lim }\limits _{N\rightarrow \infty }\left \vert \gamma _{N,0}^{(k)}-{\gamma }_{0}^{\otimes k}\right \vert _{1}^{(k)}=0\tqs \forall k\in \mathbb{N}, \end{eqnarray} \tag{ 5 }$$

that is, reduced density matrices are initially factorized. Here, γ₀ is a single spin density matrix and

$$\begin{eqnarray*} &&{\gamma }_{0}^{\otimes k}:= \mathop{\underbrace{{\gamma }_{0}\otimes {\gamma }_{0}\otimes \cdots \otimes {\gamma }_{0}{}}}\limits _{k}. \end{eqnarray*}$$

It is remarked that both the Hamiltonian and the initial state are site-symmetric, and hence the density matrix remains site-symmetric during the time evolution.

The MF theory assumes that reduced density matrices are always factorized,

$$\begin{eqnarray} &&{\gamma }_{N,t}^{(k)}\approx {\gamma }_{\mathrm{ MF},t}^{(k)}={\gamma }_{t}^{\otimes k}\quad \forall k\in \mathbb{N}. \end{eqnarray} \tag{ 6 }$$

The single spin density matrix γ_t obeys the Hartree equation,

$$\begin{eqnarray} \frac{\mathrm{d}}{\mathrm{d}t}{\gamma }_{t}&=&-\mathrm{i}{\mathcal{L}}_{1}^{(0)}{\gamma }_{t}-\mathrm{i}\hspace{0.167em} {\mathrm{Tr}}_{2}{\mathcal{L}}_{1 2}^{(V)}{\gamma }_{t}\otimes {\gamma }_{t}\nonumber\\ &=&-\mathrm{i}[{h}_{1},{\gamma }_{t}]+{\mathrm{Tr}}_{2}[{V}_{1 2},{\gamma }_{t}\otimes {\gamma }_{t}]. \end{eqnarray} \tag{ 7 }$$

The Liouville operators are defined as ${\mathcal{L}}_{i}^{(0)}(\cdot ):= [{h}_{i},(\cdot )]$ and ${\mathcal{L}}_{i j}^{(V)}(\cdot ):= [{V}_{i j},(\cdot )]$.

The statement of exactness of the MF dynamics is that if the initial state satisfies equation (5),

$$\begin{eqnarray} &&\mathop{\lim }\limits _{N\rightarrow \infty }\Vert \gamma _{N,t}^{(k)}-{\gamma }_{t}^{\otimes k}\Vert _{1}^{(k)}=0 \end{eqnarray} \tag{ 8 }$$

for any fixed k∈N and t > 0.¹

From equation (8) it is immediately verified, by using the inequality $\Vert A B\Vert _{1}^{(k)}\leq \Vert A\Vert _{\infty }\Vert B\Vert _{1}^{(k)}$, that

$$\begin{eqnarray} &&\mathop{\lim }\limits _{N\rightarrow \infty }\vert \langle {\mathcal{O}}^{(k)}\rangle _{N,t}-\langle {\mathcal{O}}^{(k)}\rangle _{\mathrm{ MF},t}\vert =0, \end{eqnarray} \tag{ 9 }$$

for any bounded operator ^(k) acting on ${\mathop{\bigotimes }\nolimits }_{i=1}^{k}{\mathcal{H}}_{i}$. Here, 〈⋅ 〉_N,t:=Tr (⋅)ρ_N,t and $\langle \cdot ~\rangle _{\mathrm{MF},t}:= {\mathrm{Tr}}_{1,2,\ldots ,k}(\cdot ){\gamma }_{\mathrm{MF},t}^{(k)}$.

We follow [22] for the proof of equation (8). We start with the Liouville equation

$$\begin{eqnarray} &&\frac{\mathrm{d}}{\mathrm{d}t}{\rho }_{N,t}=-\mathrm{i}\left (\sum _{i=1}^{N}{\mathcal{L}}_{i}^{(0)}+\frac{1}{N}\sum _{i,j=1}^{N}{\mathcal{L}}_{i j}^{(V)}\right ){\rho }_{N,t}. \end{eqnarray} \tag{ 10 }$$

By tracing out over ${\mathop{\bigotimes }\nolimits }_{i=k+1}^{N}{\mathcal{H}}_{k}$, we obtain the hierarchical equations (BBGKY hierarchy) for $\{ {\gamma }_{N,t}^{(k)}\}$ [22],

$$\begin{eqnarray} &&\frac{\mathrm{d}}{\mathrm{d}t}{\gamma }_{N,t}^{(k)}=-\mathrm{i}\sum _{i=1}^{k}{\mathcal{L}}_{i}^{(0)}{\gamma }_{N,t}^{(k)}+{\mathcal{W}}^{(k)}{\gamma }_{N,t}^{(k+1)}+\frac{1}{N}{\mathcal{V}}^{(k)}{\gamma }_{N,t}^{(k)}-\frac{k}{N}{\mathcal{W}}^{(k)}{\gamma }_{N,t}^{(k+1)}. \end{eqnarray} \tag{ 11 }$$

The super-operators ^(k) and ^(k) are defined by

$$\begin{eqnarray} &&\begin{array}{@{}l@{}} \displaystyle {\mathcal{V}}^{(k)}{\gamma }_{N,t}^{(k)}:= -\mathrm{i}\sum _{i,j=1}^{k}{\mathcal{L}}_{i j}^{(V)}{\gamma }_{N,t}^{(k)},\\ \displaystyle {\mathcal{W}}^{(k)}{\gamma }_{N,t}^{(k+1)}:= -\mathrm{i}\hspace{0.167em} {\mathrm{Tr}}_{k+1}\sum _{i=1}^{k}{\mathcal{L}}_{i,k+1}^{(V)}{\gamma }_{N,t}^{(k+1)}. \end{array} \end{eqnarray} \tag{ 12 }$$

The MF solution ${\gamma }_{\mathrm{MF},t}^{(k)}={\gamma }_{t}^{\otimes k}$ satisfies the equation

$$\begin{eqnarray} &&\frac{\mathrm{d}}{\mathrm{d}t}{\gamma }_{\mathrm{MF},t}^{(k)}=-\mathrm{i}\sum _{i=1}^{k}{\mathcal{L}}_{i}^{(0)}{\gamma }_{\mathrm{ MF},t}^{(k)}+{\mathcal{W}}^{(k)}{\gamma }_{\mathrm{ MF},t}^{(k+1)}, \end{eqnarray} \tag{ 13 }$$

which is obtained by formally taking the limit of N → ∞ in equation (11).

It should be remarked that the fact that equation (13) is obtained from equation (11) by formally taking the limit of N → ∞ alone does not ensure that the MF dynamics is exact in this limit. Because equation (11) forms a coupled chain of N equations of motion, the terms proportional to 1/N in equation (11) might be amplified in the chain and have nonnegligible contribution. In particular, for a very large k ≲ N, the last two terms on the RHS of equation (11) are no longer small. We must prove that these terms actually do not influence the dynamics of ${\gamma }_{N,t}^{(k)}$.

The next step is to construct the Duhamel series [22] by formally integrating equation (11). We define

$$\begin{eqnarray*} &&{U}_{t}^{(k)}:= \exp \left [-\mathrm{i}\sum _{i=1}^{k}{\mathcal{L}}_{i}^{(0)}t\right ]. \end{eqnarray*}$$

Then we obtain

$$\begin{eqnarray} &&\fl {\gamma }_{N,t}={U}_{t}^{(k)}{\gamma }_{N,0}^{(k)}+\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\hspace{0.167em} {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{\gamma }_{N,{t}_{1}}^{(k+1)} \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&+~\frac{1}{N}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\hspace{0.167em} {U}_{t-{t}_{1}}^{(k)}{\mathcal{V}}^{(k)}{\gamma }_{N,{t}_{1}}^{(k)}-\frac{k}{N}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\hspace{0.167em} {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{\gamma }_{N,{t}_{1}}^{(k+1)}. \end{eqnarray} \tag{ 15 }$$

We substitute this expression iteratively into the terms on the RHS not proportional to 1/N (the second term of the RHS in the above equation), and repeat this procedure L times, then we obtain the following so called Duhamel series expansion:

$$\begin{eqnarray} &&\fl {\gamma }_{N,t}^{(k)}=\sum _{l=0}^{L-1}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {\mathcal{W}}^{k+l-1}{U}_{{t}_{l}}^{(k+l)}{\gamma }_{N,0}^{(k+l)}\nonumber\\ &&+~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\hspace{0.167em} {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{L-1}-{t}_{L}}^{(k+L-1)}{\mathcal{W}}^{(k+L-1)}{\gamma }_{N,{t}_{L}}^{(k+L)}\nonumber\\ &&+~\frac{1}{N}\sum _{l=1}^{L}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{l-1}-{t}_{l}}^{(k+l-1)}{\mathcal{V}}^{(k+l-1)}{\gamma }_{N,{t}_{l}}^{(k+l-1)}\nonumber\\ &&-~\sum _{l=1}^{L}\frac{k+l-1}{N}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{l-1}-{t}_{l}}^{(k+l-1)}{\mathcal{W}}^{(k+l-1)}{\gamma }_{N,{t}_{l}}^{(k+l)}. \end{eqnarray} \tag{ 16 }$$

The Duhamel series of the MF density matrix is written as

$$\begin{eqnarray} &&\fl {\gamma }_{\mathrm{MF},t}^{(k)}=\sum _{l=0}^{L-1}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {\mathcal{W}}^{(k+l-1)}{U}_{{t}_{l}}^{(k+l)}{\gamma }_{\mathrm{ MF},0}^{(k+l)}\nonumber\\ &&+~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\hspace{0.167em} {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{L-1}-{t}_{L}}^{(k+L-1)}{\mathcal{W}}^{(k+L-1)}{\gamma }_{\mathrm{ MF},{t}_{L}}^{(k+L)}. \end{eqnarray} \tag{ 17 }$$

From these expressions, we shall show $\Vert {\gamma }_{N,t}^{(k)}-{\gamma }_{\mathrm{MF},t}^{(k)}\Vert _{1}^{(k)}\rightarrow 0$ in the limit of N → ∞. The key point is that we choose L so that 1 ≪ L ≪ N, which is realized by taking the limit of N → ∞ first and then the limit of L → ∞. From equations (16) and (17), we obtain

$$\begin{eqnarray} &&\Vert {\gamma }_{N,t}^{(k)}-{\gamma }_{\mathrm{ MF}}^{(k)}\Vert _{1}^{(k)}\leq {A}_{1}+{A}_{2}+{A}_{3}+{A}_{4}, \end{eqnarray} \tag{ 18 }$$

where the A_i (i = 1–4) are given by

$$\begin{eqnarray*} &&\fl {A}_{1}=\sum _{l=0}^{L-1}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\Vert {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {\mathcal{W}}^{k+l-1}{U}_{{t}_{l}}^{(k+l)}({\gamma }_{N,0}^{(k+l)}-{\gamma }_{\mathrm{ MF},0}^{(k+l)})\Vert _{1}^{(k)},\nonumber\\ &&\fl {A}_{2}=\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\Vert {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{L-1}-{t}_{L}}^{(k+L-1)}{\mathcal{W}}^{(k+L-1)}({\gamma }_{N,{t}_{L}}^{(k+L)}-{\gamma }_{\mathrm{ MF},{t}_{L}}^{(k+L)})\Vert _{1}^{(k)},\nonumber\\ &&\fl {A}_{3}=\frac{1}{N}\sum _{l=1}^{L}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\Vert {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{l-1}-{t}_{l}}^{(k+l-1)}{\mathcal{V}}^{(k+l-1)}{\gamma }_{N,{t}_{l}}^{(k+l-1)}\Vert _{1}^{(k)}, \end{eqnarray*}$$

and

$$\begin{eqnarray*} &&\fl {A}_{4}=\sum _{l=1}^{L}\frac{k+l-1}{N}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\Vert {U}_{t-{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1}-{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{l-1}-{t}_{l}}^{(k+l-1)}{\mathcal{W}}^{(k+l-1)}{\gamma }_{N,{t}_{l}}^{(k+l)}\Vert _{1}^{(k)}. \end{eqnarray*}$$

Now we evaluate the A_i. Since ${U}_{t}^{(k)}$ is a unitary operator, it does not change the norm, $\Vert {U}_{t}^{(k)}(\cdot )\Vert _{1}^{(k)}=\Vert (\cdot )\Vert _{1}^{(k)}$. We further use the inequality $\Vert A B\Vert _{1}^{(k)}\leq \Vert A\Vert _{\infty }\Vert B\Vert _{1}^{(k)}$. As a result, A₁ is evaluated, by recalling the definitions of ^(k) and ^(k), as

$$\begin{eqnarray*} &&{A}_{1}\leq \sum _{l=0}^{L-1}\frac{{t}^{l}}{l !}(2\Vert V\Vert _{\infty })^{l}\frac{(k+l-1)!}{(k-1)!}\left \Vert {\gamma }_{N,0}^{(k+l)}-{\gamma }_{\mathrm{ MF},0}^{(k+l)}\right \Vert _{1}^{(k+l)}. \end{eqnarray*}$$

By the assumption for the initial condition, equation (5), A₁ → 0 as N → ∞ at a fixed L. Similarly, the upper bound of A₂ becomes

$$\begin{eqnarray*} {A}_{2}&\leq &\frac{{t}^{L}}{L !}(2\Vert V\Vert _{\infty })^{L}\frac{(k+L-1)!}{(k-1)!}\left \Vert {\gamma }_{N,{t}_{L}}^{(k+l)}-{\gamma }_{\mathrm{ MF},{t}_{L}}^{(k+l)}\right \Vert _{1}^{(k+l)}\nonumber\\ &\leq &\frac{(k+L-1)!}{L !(k-1)!}(2\Vert V\Vert _{\infty }t)^{L}\times 2. \end{eqnarray*}$$

In the last inequality, we used $\left \Vert {\gamma }_{N,{t}_{L}}^{(k+l)}-{\gamma }_{\mathrm{MF},{t}_{L}}^{(k+l)}\right \Vert _{1}^{(k+l)}\leq 2$. Since (k + L − 1)!/[L!(k − 1)!] ≤ 2^k+L−1, we have

$$\begin{eqnarray*} &&{A}_{2}\leq {2}^{k}(4\Vert V\Vert _{\infty }t)^{L}. \end{eqnarray*}$$

If we restrict the time to t ≤ t₀:=1/(8||V||_∞) (this restriction will be removed later), we have A₂ ≤ 2^k−L. If we take the limit of L → ∞ after N → ∞, both A₁ and A₂ converge to zero.

Similarly, we can show that A₃ and A₄ also converge to zero in the same limit. The upper bound of A₃ is evaluated as

$$\begin{eqnarray} {A}_{3}&\leq &\frac{1}{N}\sum _{l=1}^{L-1}\frac{{t}^{l}}{l !}(2\Vert V\Vert _{\infty })^{l}(k+l-1)\frac{(k+l-1)!}{(k-1)!}\nonumber\\ &=&\frac{1}{N}\sum _{l=1}^{L-1}(2\Vert V\Vert _{\infty }t)^{l}(k+l-1)\frac{(k+l-1)!}{l !(k-1)!}\nonumber\\ &\leq &\frac{1}{N}\sum _{l=1}^{L-1}{2}^{k-1}(k+l-1)(4\Vert V\Vert _{\infty }t)^{l}. \end{eqnarray} \tag{ 19 }$$

If we take the limit of N → ∞ first, obviously A₃ goes to zero. A₄ is also bounded from above by equation (19). Therefore, it has been proven that lim_L→∞lim_N→∞A_i = 0 for the time interval 0 ≤ t ≤ t₀.

The restriction of t₀ is not essential. Since $\{ {\gamma }_{N,{t}_{0}}^{(k)}\}$ satisfies the assumption of equation (5), we can start with $\{ {\gamma }_{N,{t}_{0}}^{(k)}\}$ at time t₀ as a new initial condition. By applying the same argument as above, it is concluded that equation (8) holds in the time interval 0 ≤ t ≤ 2t₀. Repeating this argument, we can extend t₀ → ∞, and equation (8) is proved for an arbitrary time t ≥ 0.

The above proof relies on the boundedness of the operators. In optical systems, however, the above proof cannot be used as it is since the creation and annihilation operators of bosons are unbounded. We construct the proof for such a situation in the following sections.

We consider the Schrödinger dynamics of the generalized Dicke model given by

$$\begin{eqnarray} &&H(t)=\sum _{j=0}^{N}\left ({H}_{{\mathrm{ S}}_{j}}(t)+{H}_{{\mathrm{B}}_{j}}+{H}_{{\mathrm{I}}_{j}}\right )+V, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} &&\begin{array}{@{}l@{}} \displaystyle {H}_{{\mathrm{S}}_{0}}(t)={\omega }_{p}{a}^{\dagger }a+\sqrt{N}\xi (t)(a+{a}^{\dagger }),\\ \displaystyle {H}_{{\mathrm{ S}}_{j}}={\omega }_{a}\sum _{i=1}^{N}{S}_{j}^{z}\qquad (j=1,2,\ldots ,N),\\ \displaystyle {H}_{{\mathrm{ B}}_{0}}=\mathop{\sum }\limits _{r}{\omega }_{r}b_{r}^{\dagger }{b}_{r},\\ \displaystyle {H}_{{\mathrm{I}}_{0}}=\mathop{\sum }\limits _{r}{\lambda }_{r}({b}_{r}{a}^{\dagger }+b_{r}^{\dagger }a),\\ \displaystyle V=\frac{g}{\sqrt{N}}\sum _{j=1}^{N}(a{X}_{j}^{\dagger }+{a}^{\dagger }{X}_{j}). \end{array} \end{eqnarray} \tag{ 21 }$$

In H_Sj(t), the index j = 0 corresponds to the single mode of cavity photons driven by the external field ξ(t), and j = 1,2,...,N corresponds to the ensemble of N two-level atoms. Cavity photons (H_S0) interact with the free Bose field H_B0 through the coupling H_I0. Each two-level atom H_Sj, j = 1,2,...N, interacts with an arbitrary Hamilton system H_Bj through an arbitrary interaction Hamiltonian H_Ij. We call the system of the Hamiltonian H_Bj the 'j-environment' (j = 0,1,...,N), and call the composite system described by the Hamiltonian H_Sj(t) + H_Bj + H_Ij the 'j-subsystem'. We assume that the operators H_Bj + H_Ij are identical for j = 1,2,...,N; the Hamiltonian H(t) is symmetric under the exchange of two indices i and j (i, j = 1,2,...,N).

The last term of equation (20) represents the interaction between the cavity photons and the ensemble of atoms. If we choose ${X}_{j}={S}_{j}^{x}~({S}_{j}^{-})$, the Hamiltonian ${\mathop{\sum }\nolimits }_{j=0}^{N}{H}_{{\mathrm{S}}_{j}}(t)+V$ is called the Dicke model (Tavis–Cummings model [27]), with a driving force ξ. Therefore, we can regard equation (20) as an extension of the Dicke-like model; environmental systems are attached to it.

If the j-environment has a continuous spectrum, then this environment acts as a thermal reservoir in contact with H_Sj. We can obtain the Hamiltonian with a continuous spectrum as a limiting case of the discrete spectrum. It is necessary to properly choose the coupling constants {λ_r} between cavity photons and the 0-environment in order to have a well-defined limit. We introduce a parameter Λ so that the number of eigenmodes with frequencies between ω and ω + dω is given by ΛD(ω) dω. Thus the limit of Λ → ∞ corresponds to the limit of a continuum of 0-environmental modes. We assume that ${\lambda }_{r}\sim \lambda ({\omega }_{r})/\sqrt{\Lambda }.$ In this case, if we consider the limit of the continuous spectrum of the 0-environment, ${\sum }_{r}\rightarrow \Lambda \int \nolimits _{0}^{\infty }\mathrm{d}\omega \hspace{0.167em} D(\omega )$, and therefore ${\sum }_{r}\vert {\lambda }_{r}\vert ^{2}\rightarrow \int \nolimits _{0}^{\infty }\mathrm{d}\omega \hspace{0.167em} D(\omega )\lambda (\omega )^{2}=: \int \nolimits _{0}^{\infty }\mathrm{d}\omega \hspace{0.167em} J(\omega )$ does not depend on Λ. Here, J(ω) is called the spectral density [28]. We assume that

$$\begin{eqnarray} &&\mathop{\sum }\limits _{r}\vert {\lambda }_{r}\vert ^{2}\lt +\infty \end{eqnarray} \tag{ 22 }$$

in order not to diverge the effect of the 0-environment on the dynamics of cavity photons. In general, for the proof of the justification of the MF treatment, Λ is arbitrary (it can even depend on N). When Λ ≪ N, the effect of the 0-environment on the dynamics of the whole system is negligible. When Λ ∼ N, the 0-environment affects the dynamics of the total system. When Λ ≫ N, the 0-environment behaves as a thermal reservoir, and causes the dissipation.

Mean-field dynamics

The density matrix of the whole system, ρ_N,t, is defined on the Hilbert space ${\mathop{\bigotimes }\nolimits }_{i=0}^{N}{\mathcal{H}}_{i}$, where _j (j = 0,1,2,...,N) is the Hilbert space of the j-subsystem. The MF theory assumes that the k-marginal reduced density matrix ${\gamma }_{N,t}^{(k)}={\mathrm{Tr}}_{k+1,\ldots ,N}{\rho }_{\mathrm{MF},t}$ is in the product form,

$$\begin{eqnarray} &&{\gamma }_{N,t}^{(k)}\approx {\gamma }_{\mathrm{ MF},t}^{(k)}=({\rho }_{p})_{t}\otimes ({\rho }_{a})_{t}^{\otimes k}. \end{eqnarray} \tag{ 23 }$$

(ρ_p)_t is the density matrix for the 0-subsystem. Similarly, (ρ_a)_t is a common density matrix for the k-subsystem (k = 1,2,...,N). The density matrices (ρ_p)_t and (ρ_a)_t obey the following Hartree equations:

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}}{\mathrm{d}t}({\rho }_{p})_{t}=-\mathrm{i}\left [{\omega }_{p}{a}^{\dagger }a+\sqrt{N}\xi (t)(a+{a}^{\dagger })+{H}_{{\mathrm{ B}}_{0}}+{H}_{{\mathrm{I}}_{0}}\right .\nonumber\\ &&\left .+~\sqrt{N}g\left \{ [\mathrm{Tr}\hspace{0.167em} {X}_{1}({\rho }_{a})_{t}]{a}^{\dagger }+\sqrt{N}g[\mathrm{Tr}\hspace{0.167em} {X}_{1}^{\dagger }({\rho }_{a})_{t}]a\right \} ,({\rho }_{p})_{t}\right ], \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}}{\mathrm{d}t}({\rho }_{a})_{t}=-\mathrm{i}\left [{\omega }_{a}{S}_{1}^{z}+{H}_{{\mathrm{ B}}_{1}}+{H}_{{\mathrm{I}}_{1}}+\frac{g}{\sqrt{N}}\left \{ [\mathrm{Tr}\hspace{0.167em} a({\rho }_{p})_{t}]{X}_{1}^{\dagger }+[\mathrm{Tr}\hspace{0.167em} {a}^{\dagger }({\rho }_{p})_{t}]{X}_{1}\right \} ,({\rho }_{a})_{t}\right ]. \end{eqnarray} \tag{ 25 }$$

The apparent N-dependence can be removed by an appropriate scaling in the coherent state representation, see section 4.1.

It is stressed that the MF theory here corresponds to neglecting the correlation between i- and j-subsystems. In general, we cannot neglect the correlation between the system of interest (H_Sj) and the environments (H_Ij).

Rich nonequilibrium phase transitions of the model (21) have been studied with the help of the MF theory, including the optical bistability [17, 18] and the spontaneously symmetry-broken phases [19] in the regime of strong coupling (large g) and strong field (large ξ(t) = ξcos(Ωt)). Thus it is physically important to establish the validity of the MF theory in this model.

Although we cannot justify the MF theory straightforwardly in the same way as we did in the spin systems, we follow essentially the same course. We derive the BBGKY hierarchy for the reduced density matrices, construct the Duhamel series expansion, and evaluate each term. However, because of the unboundedness of boson operators we cannot use the trace norm, and the statement itself must be modified. Some preliminaries are necessary before presenting the main result of this work.

4.1. Duhamel series expansion in the coherent state representation

Since the density matrix ρ_N,t obeys the Liouville equation

$$\begin{eqnarray} &&\frac{\mathrm{d}}{\mathrm{d}t}{\rho }_{N,t}=-\mathrm{i}[H(t),{\rho }_{N,t}]=: -\mathrm{i}\mathcal{L}{\rho }_{N,t}, \end{eqnarray} \tag{ 26 }$$

we obtain the following chain of equations of motion for $\left \{ {\gamma }_{N,t}^{(k)}\right \}$:

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}}{\mathrm{d}t}{\gamma }_{N,t}^{(k)}=-\mathrm{i}\left [\sum _{j=0}^{k}\left ({H}_{{\mathrm{ S}}_{j}}+{H}_{{\mathrm{B}}_{j}}+{H}_{{\mathrm{I}}_{j}}\right ),{\gamma }_{N,t}^{(k)}\right ]\nonumber\\ &&-~\mathrm{i}\left [\frac{g}{\sqrt{N}}\left (a\sum _{j=1}^{k}{X}_{j}^{\dagger }+{a}^{\dagger }\sum _{j=1}^{k}{X}_{j}\right ),{\gamma }_{N,t}^{(k)}\right ]\nonumber\\ &&-~\mathrm{i}(N-k){\mathrm{Tr}}_{k+1}\left [\frac{g}{\sqrt{N}}\left (a{X}_{k+1}^{\dagger }+{a}^{\dagger }{X}_{k+1}\right ),{\gamma }_{N,t}^{(k+1)}\right ]. \end{eqnarray} \tag{ 27 }$$

Because it is hard to see which terms are important in the limit of N → ∞, and because the relation to the MF theory is not obvious in this form, let us introduce the Husimi 'Q-representation' of the reduced density matrix [29],

$$\begin{eqnarray} &&{Q}_{N,t}^{(k)}(\alpha ,\{ {\beta }_{r}\} )=\left (\frac{N}{\pi }\right )^{\mathcal{N}+1}{\text{ }}_{0}\langle \sqrt{N}\alpha ;\{ \sqrt{N}{\beta }_{r}\} \vert {\gamma }_{N,t}^{(k)}\vert \sqrt{N}\alpha ;\{ \sqrt{N}{\beta }_{r}\} \rangle _{0}. \end{eqnarray} \tag{ 28 }$$

Here, is the number of modes of the 0-environment ($\mathcal{N}={\sum }_{r}1\simeq \Lambda \int \nolimits _{0}^{\infty }D(\omega )\hspace{0.167em} \mathrm{d}\omega$). The vector $\vert \sqrt{N}\alpha ,\{ \sqrt{N}{\beta }_{r}\} \rangle _{0}$ denotes the coherent state on the Hilbert space ₀, and it satisfies

$$\begin{eqnarray} &&a\vert \sqrt{N}\alpha ,\{ \sqrt{N}{\beta }_{r}\} \rangle _{0}=\sqrt{N}\alpha \vert \sqrt{N}\alpha ,\{ \sqrt{N}{\beta }_{r}\} \rangle _{0}, \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&{b}_{r}\vert \sqrt{N}\alpha ,\{ \sqrt{N}{\beta }_{r}\} \rangle _{0}=\sqrt{N}{\beta }_{r}\vert \sqrt{N}\alpha ,\{ \sqrt{N}{\beta }_{r}\} \rangle _{0}. \end{eqnarray} \tag{ 30 }$$

Because ${\gamma }_{N,t}^{(k)}$ is a matrix on ${\mathop{\bigotimes }\nolimits }_{j=0}^{k}{\mathcal{H}}_{j}$ and $\left \vert \sqrt{N}\alpha ,\{ \sqrt{N}{\beta }_{r}\} \right \rangle _{0}$ is a vector on ₀, ${Q}_{N,t}^{(k)}$ is still a matrix on ${\mathop{\bigotimes }\nolimits }_{j=1}^{k}{\mathcal{H}}_{j}$, that is, the Hilbert space of k atoms and their environments. The average in equation (28) is taken only over the Hilbert space of the 0-subsystem, ₀.

${Q}_{N,t}^{(k)}$ has the following properties:

$$\begin{eqnarray*} &&\fl \text{ (i) }\quad {Q}_{N,t}^{(k)}\geq 0,\\ &&\fl \text{ (ii) }\quad \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{Q}_{N,t}^{(k)}=1,\\ &&\fl \text{ (iii) }\quad {\mathrm{Tr}}_{0,1,\ldots ,k}{\mathcal{O}}^{(k)}\left (\frac{a}{\sqrt{N}},\frac{{a}^{\dagger }}{\sqrt{N}},\left \{ \frac{{b}_{r}}{\sqrt{N}},\frac{{b}_{r}^{\dagger }}{\sqrt{N}}\right \} \right ){\gamma }_{N,t}^{(k)}\\ &&=~\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}(\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} ){Q}_{N,t}^{(k)}, \end{eqnarray*}$$

where ^(k) is an arbitrary anti-normal ordered operator acting on ${\mathop{\bigotimes }\nolimits }_{j=0}^{k}{\mathcal{H}}_{j}$; namely, all the annihilation operators are in the left of all the creation operators. From the above properties, we can obtain the expectation value of any observable ^(k) from ${Q}_{N,t}^{(k)}$.

An advantage to using ${Q}_{N,t}^{(k)}$ instead of ${\gamma }_{N,t}^{(k)}$ is that it becomes easy to see which terms are important and which terms are likely to be negligible in the limit of N → ∞. Indeed, the time evolution equation of ${Q}_{N,t}^{(k)}$ is given by

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}}{\mathrm{d}t}{Q}_{N,t}^{(k)}=-\mathrm{i}({\mathcal{L}}_{p}(t)+{\mathcal{L}}_{a}^{(k)}(\alpha )){Q}_{N,t}^{(k)}+\left (1-\frac{k}{N}\right ){\mathcal{W}}^{(k)}{Q}_{N,t}^{(k+1)}+\frac{1}{N}{\mathcal{V}}^{(k)}{Q}_{N,t}^{(k)}, \end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray} &&\fl -\mathrm{i}{\mathcal{L}}_{p}(t){Q}_{N,t}^{(k)}:= \mathrm{i}\frac{\partial }{\partial \alpha }\left [\left ({\omega }_{p}\alpha +\xi (t)+\mathop{\sum }\limits _{r}{\lambda }_{r}{\beta }_{r}\right ){Q}_{N,t}^{(k)}\right ]\nonumber\\ &&-~\mathrm{i}\frac{\partial }{\partial {\alpha }^{\ast }}\left [\left ({\omega }_{p}{\alpha }^{\ast }+\xi (t)+\mathop{\sum }\limits _{r}\lambda _{r}^{\ast }{\beta }_{r}^{\ast }\right ){Q}_{N,t}^{(k)}\right ]\nonumber\\ &&+~\mathrm{i}\mathop{\sum }\limits _{r}\left \{ \frac{\partial }{\partial {\beta }_{r}}\left [\left ({\omega }_{r}{\beta }_{r}+{\lambda }_{r}\alpha \right )Q_{N,t}^{(k)}\right ]-~\frac{\partial }{\partial \beta _{r}^{\ast }}\left [\left ({\omega }_{r}\beta _{r}^{\ast }+{\lambda }_{r}^{\ast }{\alpha }^{\ast }\right )Q_{N,t}^{(k)}\right ]\right \} , \end{eqnarray} \tag{ 32 }$$

which corresponds to the free time evolution of the 0-subsystem, and

$$\begin{eqnarray} &&\fl -\mathrm{i}{\mathcal{L}}_{a}^{(k)}(\alpha ){Q}_{N,t}^{(k)}:= -\mathrm{i}\left [\sum _{j=1}^{k}\left ({\omega }_{a}{S}_{j}^{z}+{H}_{{\mathrm{ B}}_{j}}+{H}_{{\mathrm{I}}_{j}}\right )+g\sum _{j=1}^{k}\left (\alpha {X}_{j}^{\dagger }+{\alpha }^{\ast }{X}_{j}\right ),{Q}_{N,t}^{(k)}\right ], \end{eqnarray} \tag{ 33 }$$

which represents the time evolution for k atoms and environments under the 'effective field' $g{\mathop{\sum }\nolimits }_{j=1}^{k}(\alpha {X}_{j}^{\dagger }+{\alpha }^{\ast }{X}_{j})$. The super-operators ^(k) and ^(k) are defined by

$$\begin{eqnarray} &&{\mathcal{V}}^{(k)}{Q}_{N,t}^{(k)}:= \mathrm{i}g\frac{\partial }{\partial \alpha }{Q}_{N,t}^{(k)}\sum _{j=1}^{k}{X}_{j}-\mathrm{i}g\sum _{j=1}^{k}{X}_{j}^{\dagger }\frac{\partial }{\partial {\alpha }^{\ast }}{Q}_{N,t}^{(k)}, \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} &&{\mathcal{W}}^{(k)}{Q}_{N,t}^{(k+1)}:= \mathrm{i}g\frac{\partial }{\partial \alpha }{\mathrm{Tr}}_{k+1}{X}_{k+1}{Q}_{N,t}^{(k+1)}-\mathrm{i}g\frac{\partial }{\partial {\alpha }^{\ast }}{\mathrm{Tr}}_{k+1}{X}_{k+1}^{\dagger }{Q}_{N,t}^{(k+1)}. \end{eqnarray} \tag{ 35 }$$

As was mentioned in section 3, it is remarked that if we consider the case in which the 0-environment acts as a thermal bath, the limit of Λ → ∞ should be taken before N → ∞. By assumption, {λ_r} satisfies ∑_r|λ_r|² < ∞ in the limit of Λ → ∞. In this case, the state of the 0-environment is almost unchanged, that is, $\langle {b}_{r}\rangle \sim \sqrt{N}{\beta }_{r}\sim \mathrm{O}(\sqrt{N}{\lambda }_{r})\sim \mathrm{O}(\sqrt{N/\Lambda })$, if the 0-environment is initially in equilibrium, so it is consistent with the interpretation that the 0-environment is a thermal reservoir. On the other hand, as is seen in the above estimation, when Λ ∼ O(N), the state of the 0-environment is strongly disturbed, and it cannot be regarded as a thermal bath in this case. Anyway, the justification of the MF approach is possible for both these two cases.

If we formally take the limit of N → ∞ in equation (31), we obtain the following equation of motion:

$$\begin{eqnarray} &&\frac{\mathrm{d}}{\mathrm{d}t}{Q}_{\mathrm{MF},t}^{(k)}=-\mathrm{i}({\mathcal{L}}_{p}(t)+{\mathcal{L}}_{a}^{(k)}(\alpha )){Q}_{\mathrm{ MF},t}^{(k)}+{\mathcal{W}}^{(k)}{Q}_{\mathrm{ MF},t}^{(k+1)}. \end{eqnarray} \tag{ 36 }$$

If initially the density matrix is in the product form,

$$\begin{eqnarray*} &&{Q}_{\mathrm{MF},0}^{(k)}={Q}_{\mathrm{ MF},0}^{(0)}\otimes ({\rho }_{a})_{0}^{\otimes k},\tqs \forall k\in \mathbb{N}, \end{eqnarray*}$$

the solution of equation (36) is equivalent to that of the MF dynamical equation. Therefore, we shall compare the solutions of equation (31) with those of equation (36).

We define the time evolution operator due to $-\mathrm{i}{\mathcal{L}}_{\mathrm{p}}(t)-\mathrm{i}{\mathcal{L}}_{a}^{(k)}(\alpha )$ as

$$\begin{eqnarray} &&{U}_{t,s}^{(k)}:= \overleftarrow {\mathrm{T}}\exp \left [-\mathrm{i}\int \nolimits _{s}^{t}\mathrm{d}{t}^{\prime}\left ({\mathcal{L}}_{p}({t}^{\prime})+{\mathcal{L}}_{a}^{(k)}(\alpha )\right )\right ]. \end{eqnarray} \tag{ 37 }$$

Here $\overleftarrow {\mathrm{T}}$ is the time-ordering operator (the arrow implies the direction from past to future). Then we can obtain the formal integral equation for ${Q}_{N,t}^{(k)}$ from equation (31) as

$$\begin{eqnarray} &&\fl {Q}_{N,t}^{(k)}={U}_{t,0}^{(k)}{Q}_{N,0}^{(k)}(0)+\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{Q}_{N,{t}_{1}}^{(k+1)}\nonumber\\ &&+~\frac{1}{N}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{V}}^{(k)}{Q}_{N,{t}_{1}}^{(k)}-\frac{k}{N}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{Q}_{N,{t}_{1}}^{(k+1)}. \end{eqnarray} \tag{ 38 }$$

It is almost the same as equation (11). By substituting this expression into the RHS iteratively L times, we obtain the following Duhamel series expansion:

$$\begin{eqnarray} &&\fl {Q}_{N,t}^{(k)}=\sum _{l=0}^{L-1}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {\mathcal{W}}^{k+l-1}{U}_{{t}_{l},0}^{(k+l)}{Q}_{N,0}^{(k+l)}\nonumber\\ &&+~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{L-1},{t}_{L}}^{(k+L-1)}{\mathcal{W}}^{(k+L-1)}{Q}_{N,{t}_{L}}^{(k+L)}\nonumber\\ &&+~\frac{1}{N}\sum _{l=1}^{L}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}{\mathcal{V}}^{(k+l-1)}{Q}_{N,{t}_{l}}^{(k+l-1)}\nonumber\\ &&-~\sum _{l=1}^{L}\frac{k+l-1}{N}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}{\mathcal{W}}^{(k+l-1)}{Q}_{N,{t}_{l}}^{(k+l)}. \end{eqnarray} \tag{ 39 }$$

Similarly, we can obtain the Duhamel expansion for the MF density matrix ${Q}_{\mathrm{MF},t}^{(k)}$,

$$\begin{eqnarray} &&\fl {Q}_{\mathrm{MF},t}^{(k)}=\sum _{l=0}^{L-1}\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {\mathcal{W}}^{k+l-1}{U}_{{t}_{l},0}^{(k+l)}{Q}_{\mathrm{ MF},0}^{(k+l)}\nonumber\\ &&+~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\cdots {U}_{{t}_{L-1},{t}_{L}}^{(k+L-1)}{\mathcal{W}}^{(k+L-1)}{Q}_{\mathrm{ MF},{t}_{L}}^{(k+L)}. \end{eqnarray} \tag{ 40 }$$

These expansions are the starting point for justification of the MF theory.

4.2. Restriction of observables

Let ^(k) be an operator acting on the Hilbert space ${\mathop{\bigotimes }\nolimits }_{j=0}^{k}{\mathcal{H}}_{j}$. Because the creation and annihilation operators for bosons are unbounded operators, it is hard to prove that the expectation value of an arbitrary operator ^(k) calculated by ${Q}_{N,t}^{(k)}$ coincides with that calculated by Q_MF,t. Therefore, we now restrict the physical quantities.

In this paper, we focus only on the observables ${\mathcal{O}}^{(k)}(\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} )\in \mathcal{B}$. The set of observables ℬ is defined as follows. All the observables ${\mathcal{O}}^{(k)}(\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} )\in \mathcal{B}$ satisfy the following conditions. There exist some complex numbers c₀,c_a,{c_r}∈C, with lim_Λ→∞∑_r|c_r|² <+ ∞, and some positive integer s > 0 and some positive numbers κ > 0 and {a_q > 0}, q = 0,1,...,s, such that

$$\begin{eqnarray} &&\left \Vert \frac{{\partial }^{m}{\mathcal{O}}^{(k)}(\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} )}{\partial {\alpha }_{{\sigma }_{1}}\partial {\alpha }_{{\sigma }_{2}}\cdots \partial {\alpha }_{{\sigma }_{m}}}\right \Vert _{\infty }\leq {\kappa }^{m}\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha +\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}\right \vert ^{q}, \end{eqnarray} \tag{ 41 }$$

for all m = 0,1,... and for all $\{ {\sigma }_{j}\} _{j=1}^{m}$. Here, σ_j =± 1, α₁:=α, and α₋₁:=α*. For instance, all the operators expressed by a polynomial of $\{ a,{a}^{\dagger },\{ {\vec{S}}_{i}\} _{i=1}^{k}\}$ belong to ℬ after the correspondence $\{ a/\sqrt{N}\rightarrow \alpha ,{a}^{\dagger }/\sqrt{N}\rightarrow {\alpha }^{\ast }\}$ is made.

The derivatives ∂/∂α and ∂/∂α* appear because of the quantum fluctuation, i.e. the commutation relation [a,a^†] = 1. Roughly speaking, the condition (41) means that we exclude observables that are too sensitive to quantum fluctuations. For example, an operator like exp[a^†a/N], which roughly corresponds to ^(k) = exp[α*α], is excluded for any value of  > 0 although its expectation value might exist and be finite.

4.3. Free time evolution of the 0-subsystem

In this section, we introduce several time evolution operators and discuss their relation. The relation given by equation (45) corresponds to the transformation from the 'Schrödinger picture' to the 'interaction picture' in terms of the coherent state representation.

Let us define

$$\begin{eqnarray} &&{U}_{t,s}^{(p)}:= \overleftarrow {\mathrm{T}}\exp \left [-\mathrm{i}\int \nolimits _{s}^{t}\mathrm{d}{t}^{\prime}\hspace{0.167em} {\mathcal{L}}_{p}({t}^{\prime})\right ]. \end{eqnarray} \tag{ 42 }$$

For an arbitrary function $f(\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} )$,

$$\begin{eqnarray} &&f(\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} ){U}_{t,0}^{(p)}={U}_{t,0}^{(p)}f(\alpha (t),{\alpha }^{\ast }(t),\{ {\beta }_{r}(t),{\beta }_{r}^{\ast }(t)\} ), \end{eqnarray} \tag{ 43 }$$

where α(t) is the solution of the following equations:

$$\begin{eqnarray} &&\begin{array}{@{}l@{}} \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}\alpha (t)=-\mathrm{i}\left [{\omega }_{p}\alpha (t)+\xi (t)+\mathop{\sum }\limits _{r}{\lambda }_{r}{\beta }_{r}(t)\right ],\\ \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}{\beta }_{r}(t)=-\mathrm{i}\left ({\omega }_{r}{\beta }_{r}(t)+{\lambda }_{r}\alpha (t)\right ), \end{array} \end{eqnarray} \tag{ 44 }$$

with the initial conditions α(0) = α and β_r(0) = β_r.

The following property derived from equation (43) is important:

$$\begin{eqnarray} &&{U}_{t,s}^{(k)}{U}_{s,0}^{(p)}={U}_{t,0}^{(p)}{\tilde {U}}_{t,s}^{(k)}, \end{eqnarray} \tag{ 45 }$$

where

$$\begin{eqnarray} &&{\tilde {U}}_{t,s}^{(k)}:= \overleftarrow {T}\exp \left [-\mathrm{i}\int \nolimits _{s}^{t}\mathrm{d}{t}^{\prime}\hspace{0.167em} {\mathcal{L}}_{a}^{(k)}(\alpha ({t}^{\prime}))\right ]. \end{eqnarray} \tag{ 46 }$$

This relation allows us to divide the time evolution operator ${U}_{t,s}^{(k)}$ into two parts, the free time evolution of the 0-subsystem ${U}_{t,0}^{(p)}$ and the time evolution of the remaining part ${\tilde {U}}_{t,s}^{(k)}$. This is viewed as the transformation to the 'interaction picture' in the sense that the photon amplitude α(t) appearing in k atoms' time evolution operator ${\tilde {U}}_{t,s}^{(k)}$ evolves under the Hamiltonian of the 0-subsystem.

Because these equations are linear, the solutions are written in the form

$$\begin{eqnarray} &&\begin{array}{@{}l@{}} \displaystyle \alpha (t)={h}_{a}(t)+{g}_{a a}(t)\alpha +\mathop{\sum }\limits _{r}{g}_{a r}(t){\beta }_{r},\\ \displaystyle {\beta }_{r}(t)={h}_{r}(t)+{g}_{r a}(t)\alpha +\mathop{\sum }\limits _{{r}^{\prime}}{g}_{r{r}^{\prime}}(t){\beta }_{{r}^{\prime}}. \end{array} \end{eqnarray} \tag{ 47 }$$

It is noted that

$$\begin{eqnarray} &&\begin{array}{@{}l@{}} \displaystyle \vert {g}_{a a}(t)\vert ^{2}+\mathop{\sum }\limits _{r}\left \vert {g}_{a r}(t)\right \vert ^{2}=1,\\ \displaystyle \vert {g}_{r a}(t)\vert ^{2}+\mathop{\sum }\limits _{{r}^{\prime}}\left \vert {g}_{r{r}^{\prime}}(t)\right \vert ^{2}=1. \end{array} \end{eqnarray} \tag{ 48 }$$

The matrix g(t) is a unitary and symmetric matrix. This property will be used later. Here, it should be noted that it is an important assumption that the 0-subsystem is a linear Bose system.

From now on, we prove that the MF theory is exact in the limit of N → ∞ in the sense that if initially

$$\begin{eqnarray} &&\fl \mathop{\lim }\limits _{N\rightarrow \infty }\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\left ({Q}_{N,0}^{(k)}-{Q}_{\mathrm{ MF},0}^{(k)}\right )=0 \end{eqnarray} \tag{ 49 }$$

for ∀^(k)∈ℬ and any fixed k∈N, it implies that

$$\begin{eqnarray} &&\fl \mathop{\lim }\limits _{N\rightarrow \infty }\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\left ({Q}_{N,t}^{(k)}-{Q}_{\mathrm{ MF},t}^{(k)}\right )=0 \end{eqnarray} \tag{ 50 }$$

for any fixed time t > 0.

5.1. Finiteness of expectation values of observables

First of all, we prove that all the expectation values of observables ^(k)∈ℬ at time t are finite as long as they are also finite at initial time t = 0. Unboundedness of operators is obstructive for our proof, hence this property is desirable. Once we can show this property, the justification of the MF approximation is almost straightforward.

From equation (39), we obtain

$$\begin{eqnarray} &&\fl \left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right ){Q}_{N,t}^{(k)}\right \vert \leq {B}_{1}+{B}_{2}+{B}_{3}+{B}_{4}. \end{eqnarray} \tag{ 51 }$$

The {B_i} are defined by

$$\begin{eqnarray*} &&\fl {B}_{1}=\sum _{l=0}^{L-1}\left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\right .\nonumber\\ &&\fl \quad \left .\times \int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {\mathcal{W}}^{(k+l-1)}{U}_{{t}_{l},0}^{(k+l)}{Q}_{N,0}^{(k+l)}\right \vert ,\\ &&\fl {B}_{2}=\left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\right .\nonumber\\ &&\fl \quad \left .\vphantom {\int {\mathrm {d}}^{2}\alpha \left (\mathop {\prod } _{r}\int {\mathrm {d}}^{2}{\beta }_{r}\right ){\mathrm {Tr}}_{1,2,\ldots ,k}{\mathcal {O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )} \times ~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {U}_{{t}_{L-1},{t}_{L}}^{(k+L-1)}{\mathcal{W}}^{(k+L-1)}{Q}_{N,{t}_{L}}^{(k+L)}\right \vert ,\\ &&\fl {B}_{3}=\frac{1}{N}\sum _{l=1}^{L}\left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\right .\nonumber\\ &&\fl \quad \left .\vphantom {\vert \int {\mathrm {d}}^{2}\alpha \left (\mathop {\prod } _{r}\int {\mathrm {d}}^{2}{\beta }_{r}\right ){\mathrm {Tr}}_{1,2,\ldots ,k}{\mathcal {O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )} \times ~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {U}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}{\mathcal{V}}^{(k+l-1)}{Q}_{N,{t}_{l}}^{(k+l-1)}\right \vert , \end{eqnarray*}$$

and

$$\begin{eqnarray*} &&\fl {B}_{4}=\sum _{l=1}^{L}\frac{k+l-1}{N}\left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\right .\nonumber\\ &&\fl \quad \left .\vphantom {\vert \int {\mathrm {d}}^{2}\alpha \left (\mathop {\prod } _{r}\int {\mathrm {d}}^{2}{\beta }_{r}\right ){\mathrm {Tr}}_{1,2,\ldots ,k}{\mathcal {O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )} \times ~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {U}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}{\mathcal{W}}^{(k+l-1)}{Q}_{N,{t}_{l}}^{(k+l)}\right \vert . \end{eqnarray*}$$

What we have to do is to evaluate the upper bounds of B_i (i = 1,...,4). Since the derivation is very complicated, we give the derivation in appendix A, and here we just mention the strategy briefly. The analysis is similar to that in section 2, but due to the derivatives with respect to α and α* in ^(k), equation (34), and ^(k), equation (35), it must be modified. First, we move to the 'interaction picture' by using equation (45). By integrating by part repeatedly, we rewrite {B_i} so that all the derivatives are acting only on $\{ {\tilde {U}}_{{t}_{n},{t}_{n+1}}^{(k+n)}\}$ (n = 1,2,...,l or L) and ^(k). Then we can show that the derivatives of ${\tilde {U}}_{{t}_{n},{t}_{n+1}}^{(k+n)}$ with respect to α and α* are bounded above by equation (A.7). In addition, the derivatives of ^(k) are also bounded due to the restriction of observables, see equation (41). By using these bounds, we can obtain the upper bounds of {B_i} which approach zero in the limit of 1 ≪ L ≪ N.

The derived upper bounds are the following:

$$\begin{eqnarray} &&\fl {B}_{1}\leq {2}^{k}{\mathrm{e}}^{2 \kappa }\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,0}^{(0)}, \end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray} &&\fl {B}_{2}\leq {2}^{-L+k-1}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}^{\prime}}^{(0)}, \end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray} &&\fl {B}_{3}\lt \frac{1}{N}{2}^{k-1}(k+2){\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}^{\prime}}^{(0)}.\nonumber\\ \end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray} &&\fl {B}_{4}\leq \frac{1}{N}{2}^{k-1}(k+1){\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}^{\prime}}^{(0)}.\nonumber\\ \end{eqnarray} \tag{ 55 }$$

From the derived upper bounds of B_i (i = 1–4), we show that the expectation value of any observable in ℬ at any fixed time t is finite. Collecting the derived upper bounds presented above, we obtain

$$\begin{eqnarray} &&\fl \vert \langle {\mathcal{O}}^{(k)}\rangle _{t}\vert \lt {2}^{k}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\left \langle \left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}\right \rangle _{0}\nonumber\\ &&+~{2}^{k-1}{\mathrm{e}}^{2 \kappa }\left ({2}^{-L}+\frac{2 k+3}{N}\right )\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\left \langle \left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}\right \rangle _{{t}^{\prime}}\nonumber\\ &&\equiv ~{2}^{k}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\left \langle \left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}\right \rangle _{0}\nonumber\\ &&+~\sum _{q=0}^{s}{b}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\left \langle \left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}\right \rangle _{{t}^{\prime}}, \end{eqnarray} \tag{ 56 }$$

where b_q = 2^k−1e^2κ(2^−L + (2k + 3)/N)a_q. We can see that b_q → 0 in the limit of N → ∞ and L → ∞. Therefore, in the limit of N → ∞, we can choose b_q as an arbitrarily small value. Here, $\langle \cdot ~\rangle _{t}:= \int {\mathrm{d}}^{2}\alpha \hspace{0.167em} (\prod \int {\mathrm{d}}^{2}{\beta }_{r})(\cdot ){Q}_{N,t}^{(0)}$ is shorthand notation for the expectation value. In order to show the finiteness, we must show that

$$\begin{eqnarray*} &&\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\left \langle \left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}\right \rangle _{{t}^{\prime}} \end{eqnarray*}$$

is finite.

In order to do so, we put

$$\begin{eqnarray} &&{\mathcal{O}}^{(k)}=\left \vert {d}_{0}+{d}_{a}\alpha +\mathop{\sum }\limits _{r}{d}_{r}{\beta }_{r}\right \vert ^{n}=: {\mathcal{O}}_{n}, \end{eqnarray} \tag{ 57 }$$

and define

$$\begin{eqnarray} &&{\mathcal{O}}_{n}(t):= \left \vert {d}_{0}+{d}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{d}_{r}{\beta }_{r}(t)\right \vert ^{n}. \end{eqnarray} \tag{ 58 }$$

Because the constants d₀, d_a, and ∑_r|d_r|² are finite, the generality is not lost if we restrict these constants as

$$\begin{eqnarray} &&\vert \vec{d}\vert ^{2}:= \vert {d}_{0}\vert ^{2}+\vert {d}_{a}\vert ^{2}+\mathop{\sum }\limits _{r}\vert {d}_{r}\vert ^{2}\leq 1. \end{eqnarray} \tag{ 59 }$$

Obviously, _n∈ℬ with s = n, c₀ = d₀, c_a = d_a, c_r = d_r, κ = nd_a, and a_q = 1.

Now we evaluate the quantity

$$\begin{eqnarray*} &&{\mathcal{X}}_{n}:= \mathop{\max }\limits _{\vec{d}:\vert \vec{d}\vert \leq 1}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\langle {\mathcal{O}}_{n}\rangle _{{t}^{\prime}}. \end{eqnarray*}$$

From equation (56), we have

$$\begin{eqnarray} &&{\mathcal{X}}_{n}\lt {2}^{k}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{n}\langle {\mathcal{O}}_{q}(t)\rangle _{0}+b\sum _{q=0}^{n}\mathop{\max }\limits _{\vec{d}:\vert \vec{d}\vert \leq 1}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\langle {\mathcal{O}}_{q}(t)\rangle _{{t}^{\prime}}, \end{eqnarray} \tag{ 60 }$$

where b = 2^k−1e^2κ[2^−L + (2k + 3)/N].

Now we use the properties of α(t) and β_r(t). By substituting equation (47) into the definition of _q(t), we have

$$\begin{eqnarray} &&{\mathcal{O}}_{q}(t)=\left \vert {d}_{0}^{\prime}+{d}_{a}^{\prime}\alpha +\mathop{\sum }\limits _{r}d_{r}^{\prime}{\beta }_{r}\right \vert ^{q}, \end{eqnarray} \tag{ 61 }$$

where

$$\begin{eqnarray*} &&{d}_{0}^{\prime}={d}_{0}+{d}_{a}{h}_{a}(t)+\mathop{\sum }\limits _{r}{d}_{r}{h}_{r}(t)=: \vec{d}\cdot \vec{h}(t),\\ &&d_{a}^{\prime}={d}_{a}{g}_{a a}(t)+\mathop{\sum }\limits _{r}{d}_{r}{g}_{r a}(t),\\ &&d_{r}^{\prime}={d}_{a}{g}_{a r}(t)+\mathop{\sum }\limits _{{r}^{\prime}}{d}_{{r}^{\prime}}{g}_{{r}^{\prime}r}(t). \end{eqnarray*}$$

The norm of this new vector $\vec{{d}^{\prime}}$ is given by

$$\begin{eqnarray*} &&\vert \vec{{d}^{\prime}}\vert ^{2}=\vert {d}_{0}^{\prime}\vert ^{2}+\vert {d}_{a}^{\prime}\vert ^{2}+\mathop{\sum }\limits _{r}\vert d_{r}^{\prime}\vert ^{2}\leq \vert \vec{h}(t)\vert ^{2}+\mathop{\sum }\limits _{x}\left \vert \mathop{\sum }\limits _{{x}^{\prime}}{d}_{{x}^{\prime}}{g}_{{x}^{\prime}x}(t)\right \vert ^{2}, \end{eqnarray*}$$

where x,x'∈{a,{r}}. Now we use the fact that the matrix g(t) is unitary (the absolute value of any eigenvalue of g(t) is unity). This yields

$$\begin{eqnarray} &&\vert \vec{{d}^{\prime}}\vert ^{2}\leq \vert \vec{h}(t)\vert ^{2}+1, \end{eqnarray} \tag{ 62 }$$

and the RHS is finite for any fixed time t because of the linearity of equation (44) and the assumption of non-pathological external field. This also indicates that _q(t)∈ℬ, therefore, the first term of the RHS in equation (60) is finite. We write it as ${f}_{n}(t)={2}^{k}{\mathrm{e}}^{2 \kappa }{\mathop{\sum }\nolimits }_{q=0}^{n}\left \langle {\mathcal{O}}_{q}(t)\right \rangle _{0}$.

We define the vector $\vec{{d}^{\prime\prime}}:= \vec{{d}^{\prime}}/(\vert \vec{h}(t)\vert ^{2}+1)$, whose norm is less than or equal to unity, $\vert \vec{{d}^{\prime\prime}}\vert \leq 1$ regardless of the choice of $\vec{d}$. Obviously $\left \{ \vec{{d}^{\prime\prime}}:\vert \vec{{d}^{\prime\prime}}\vert \leq 1\right \} \supseteq \left \{ \vec{{d}^{\prime\prime}}:\vert \vec{d}\vert \leq 1\right \}$. Therefore, from equations (60) and (61), we obtain

$$\begin{eqnarray} &&\fl {\mathcal{X}}_{n}\lt {f}_{n}(t)+b\left (\vert \vec{h}(t)\vert ^{2}+1\right )\sum _{q=0}^{n}\mathop{\max }\limits _{\vec{d}:\vert \vec{d}\vert \leq 1}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\left \langle \left \vert d_{0}^{\prime\prime}+d_{a}^{\prime\prime}\alpha +\mathop{\sum }\limits _{r}d_{r}^{\prime\prime}{\beta }_{r}\right \vert ^{q}\right \rangle _{{t}^{\prime}}\nonumber\\ &&\leq ~{f}_{n}(t)+b\left (\vert \vec{h}(t)\vert ^{2}+1\right )\sum _{q=0}^{n}\mathop{\max }\limits _{\vec{{d}^{\prime\prime}}:\vert \vec{{d}^{\prime\prime}}\vert \leq 1}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\left \langle \left \vert d_{0}^{\prime\prime}+d_{a}^{\prime\prime}\alpha +\mathop{\sum }\limits _{r}d_{r}^{\prime\prime}{\beta }_{r}\right \vert ^{q}\right \rangle _{{t}^{\prime}}\nonumber\\ &&={f}_{n}(t)+b\left (\vert \vec{h}(t)\vert ^{2}+1\right )\sum _{q=0}^{n}{\mathcal{X}}_{q}. \end{eqnarray} \tag{ 63 }$$

Since we can choose b as an arbitrarily small number, by starting from ₀ = 1, we can inductively show that all the _n are finite by using the above inequality. This completes the proof of the finiteness of the expectation values, because the RHS of equation (56) is then finite.

In particular, we find

$$\begin{eqnarray} &&\mathop{\lim }\limits _{L\rightarrow \infty }\mathop{\lim }\limits _{N\rightarrow \infty }{B}_{2}=\mathop{\lim }\limits _{L\rightarrow \infty }\mathop{\lim }\limits _{N\rightarrow \infty }{B}_{3}=\mathop{\lim }\limits _{L\rightarrow \infty }\mathop{\lim }\limits _{N\rightarrow \infty }{B}_{4}=0. \end{eqnarray} \tag{ 64 }$$

5.2. Justification of the mean-field approach

The proof proceeds very similarly to section 5.1. From equations (39) and (40), we obtain the upper bound

$$\begin{eqnarray} &&\fl \left \vert \int {\mathrm{d}}^{2}\alpha \hspace{0.167em} \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\left ({Q}_{N,t}^{(k)}-{Q}_{\mathrm{ MF},t}^{(k)}\right )\right \vert \nonumber\\ &&\leq ~{B}_{1}^{\prime}+{B}_{2}^{\prime}+{B}_{3}+{B}_{4}. \end{eqnarray} \tag{ 65 }$$

Here, B₃ and B₄ are the same as in section 5.1. ${B}_{1}^{\prime}$ and ${B}_{2}^{\prime}$ are defined by

$$\begin{eqnarray} &&\fl {B}_{1}^{\prime}:= \sum _{l=0}^{L-1}\left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\right .\nonumber\\ &&\times ~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\left .\cdots {\mathcal{W}}^{k+l-1}{U}_{{t}_{l},0}^{(k+l)}\left ({Q}_{N,0}^{(k+l)}-{Q}_{\mathrm{ MF},0}^{(k+l)}\right )\vphantom {\sum _{1}^{2}} \right \vert ,\nonumber\\ &&\fl {B}_{2}^{\prime}:= \left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\right .\nonumber\\ &&\times ~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\left .\cdots {U}_{{t}_{L-1},{t}_{L}}^{(k+L-1)}{\mathcal{W}}^{(k+L-1)}\left ({Q}_{N,{t}_{L}}^{(k+L)}-{Q}_{\mathrm{ MF},{t}_{L}}^{(k+L)}\right )\vphantom {\sum _{1}^{2}} \right \vert .\nonumber\\ \end{eqnarray} \tag{ 66 }$$

In appendix A, it is shown that lim_N→∞B₃ = lim_N→∞B₄ = 0 in the time interval t∈[0,t₀]. The proof of

$$\begin{eqnarray} &&\mathop{\lim }\limits _{L\rightarrow \infty }\mathop{\lim }\limits _{N\rightarrow \infty }B_{1}^{\prime}=\mathop{\lim }\limits _{L\rightarrow \infty }\mathop{\lim }\limits _{N\rightarrow \infty }B_{2}^{\prime}=0 \end{eqnarray} \tag{ 67 }$$

for t∈[0,t₀] is similar to the analysis in appendix A. We give the derivation of equation (67) in appendix B.

Up to now, we have shown that ${B}_{1}^{\prime}+{B}_{2}^{\prime}+{B}_{3}+{B}_{4}\rightarrow 0$ as N → ∞ for 0 ≤ t ≤ t₀. This means that equation (50) has been shown for 0 ≤ t ≤ t₀. Finally, for completeness of the proof, we must extend the time interval from t∈[0,t] to t∈[0,∞), but there is no difficulty on this point. If we regard t = t₀ as a new initial time and repeat the argument, we can show that equation (50) is correct for any fixed time t > 0. Thus the proof of the justification of the MF dynamics has been completed.

Finally we conclude the present work by making some remarks on our result.

• Our result is quite general: environments may be attached to the system of interest, and a time-dependent driving field may be applied. When an environment is large enough, it acts as a thermal reservoir on the system of interest, hence the result of this work is applicable to open quantum systems. However, we assumed a special Hamiltonian for the 0-environment, H_B0 + H_I0 in equation (21). It is preferable to generalize our result to a wider class of 0-environmental Hamiltonians in the future.
• We restricted the class of physical quantities to ℬ in this work. This restriction, however, might be just a technical assumption. Extension to a more general class of operators is an important issue. Another important issue is to treat the fluctuations. For instance, ${a}^{\dagger }a/N=(a{a}^{\dagger }-1)/N=(a/\sqrt{N})({a}^{\dagger }/\sqrt{N})-1/N$ belongs to ℬ, but $({a}^{\dagger }a-\langle {a}^{\dagger }\rangle \langle a\rangle )/\sqrt{N}$, whose expectation value converges to a finite value as N → ∞, does not. In the seminal work by Hepp and Lieb [16], the authors showed that, under the singular reservoir limit and the approximation of replacing the bosons by the ensemble of fermions, the equations of motion for the fluctuations are equal to the linearized equations of the intensive observables; this is known as Onsager's regression hypothesis [30]–[32]. We have not been able to confirm the validity of this hypothesis without those approximations yet.
• Without the driving force, and if the thermal reservoirs have identical temperatures, the MF dynamics usually predicts thermalization of the system of interest. This means that the limit of N → ∞ and the limit of t → ∞ are interchangeable under such an equilibrium situation. In other words, there is no quasi-stationary state in an open system under the equilibrium situation, which is in contrast to closed systems, see [33] for this aspect. However, it is less obvious whether or not the thermodynamic limit and the long-time limit are interchangeable under the nonequilibrium situation studied in this paper. In this work, we always take the limit of N → ∞ first and then the limit of t → ∞. In a real experimental setup, there can be a situation in which the thermodynamic limit ought to be taken after the long-time limit since the number of atoms in the cavity is not large enough.
• We can justify the MF theory, but it is a separate issue whether the MF dynamical equations, equations (24) and (25), can be solved. Since the degrees of freedom of the environmental systems remain, it is difficult to solve the MF dynamical equations exactly. It would be very interesting if we could exactly solve the quantum dynamics of the whole system including thermal reservoirs with the help of the MF theory. If this were done, we would be able to obtain some insight into the effects of non-Markovian dynamics [25] and the accuracy of the perturbative quantum master equation at long times [34, 35] in an interacting many body system. When the system of interest is small, which is the standard setting of open quantum systems, it is recognized that the non-Markovian effect is negligible in the van Hove limit [36]. When the system of interest is large, however, it is not obvious whether the use of the Born–Markov quantum master equation is justified even for the small coupling between the system of interest and the thermal bath, because two limits, the thermodynamic limit and the van Hove limit, are involved. Therefore, it is important to obtain the exact time evolution in a macroscopic open system. It is noted that Merkli and Berman [23] recently obtained a rigorous result in this direction for a simple model which is purely dephasing and where relaxation of energy does not take place. In our setting, if we assume that each spin interacts with infinitely many reservoirs of the same structure, the dynamics will be solved at least in a numerically exact manner, because we can also apply the MF approximation to the coupling between the system of interest and the reservoirs in that case.

We hope that the result of this work will become a good starting point to consider the above remaining problems and analyze some interesting nonequilibrium phenomena.

The author thanks Seiji Miyashita and Tatsuhiko Shirai for fruitful discussion. This work is supported by the Sumitomo Foundation (grant No. 120753).

A.1. Evaluation of B₁

First, we consider the term B₁. From equations (34), (35), and (45), we obtain the following expression:

$$\begin{eqnarray} &&\fl {B}_{1}=\sum _{l=0}^{L-1}\left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}(\alpha (t),{\alpha }^{\ast }(t))\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\right .\nonumber\\ &&\times ~\sum _{\{ {\sigma }_{j}=\pm 1\} }^{l}{\tilde {U}}_{t,{t}_{1}}^{(k)}\frac{\partial }{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}{\mathrm{Tr}}_{k+1}{X}_{k+1}^{({\sigma }_{1})}{\tilde {U}}_{{t}_{1},{t}_{2}}^{(k+1)}\frac{\partial }{\partial {\alpha }_{{\sigma }_{2}}({t}_{2})}{\mathrm{Tr}}_{k+2}{X}_{k+2}^{({\sigma }_{2})}\nonumber\\ &&\left .\cdots \frac{\partial }{\partial {\alpha }_{{\sigma }_{l}}({t}_{l})}{\mathrm{Tr}}_{k+l}{X}_{k+l}^{({\sigma }_{l})}{\tilde {U}}_{{t}_{l},0}^{(k+l)}{Q}_{N,0}^{(k+l)}\vphantom {\sum _{1}^{2}} \right \vert , \end{eqnarray} \tag{ A.1 }$$

where ${X}_{j}^{(1)}:= {X}_{j}$ and ${X}_{j}^{(-1)}:= {X}_{j}^{\dagger }$. By performing integration by part, we obtain

$$\begin{eqnarray} &&\fl {B}_{1}=\sum _{l=0}^{L-1}\Biggl \vert \int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k+l}\nonumber\\ &&\times ~\sum _{\{ {\sigma }_{j}=\pm 1\} }^{l}{X}_{k+1}^{({\sigma }_{1})}{X}_{k+2}^{({\sigma }_{2})}\cdots {X}_{k+l}^{({\sigma }_{l})}\nonumber\\ &&\times ~\frac{\partial }{\partial {\alpha }_{{\sigma }_{l}}({t}_{l})}\Biggl \{ \frac{\partial }{\partial {\alpha }_{{\sigma }_{l-1}}({t}_{l-1})}\Biggl \{ \cdots \frac{\partial }{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}\Biggl \{ {\mathcal{O}}^{(k)}(\alpha (t),{\alpha }^{\ast }(t),\{ {\beta }_{r}(t),{\beta }_{r}^{\ast }(t)\} )\nonumber\\ &&\times ~{\tilde {U}}_{t,{t}_{1}}^{(k)}\Biggr \} {\tilde {U}}_{{t}_{1},{t}_{2}}^{(k+1)}\Biggr \} \cdots {\tilde {U}}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}\Biggr \} {\tilde {U}}_{{t}_{l},0}^{(k+l)}{Q}_{N,0}^{(k+l)}\Biggr \vert . \end{eqnarray} \tag{ A.2 }$$

We put

$$\begin{eqnarray} &&\fl {B}_{1}^{(k+l)}:= \int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\sum _{\{ {\sigma }_{j}=\pm 1\} }^{l}{X}_{k+1}^{({\sigma }_{1})}{X}_{k+2}^{({\sigma }_{2})}\cdots {X}_{k+l}^{({\sigma }_{l})}\nonumber\\ &&\times ~\frac{\partial }{\partial {\alpha }_{{\sigma }_{l}}({t}_{l})}\Biggl \{ \frac{\partial }{\partial {\alpha }_{{\sigma }_{l-1}}({t}_{l-1})}\Biggl \{ \cdots \frac{\partial }{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}\Biggl \{ \nonumber\\ &&\times ~{\mathcal{O}}^{(k)}(\alpha (t),{\alpha }^{\ast }(t),\{ {\beta }_{r}(t),{\beta }_{r}^{\ast }(t)\} ){\tilde {U}}_{t,{t}_{1}}^{(k)}\Biggr \} {\tilde {U}}_{{t}_{1},{t}_{2}}^{(k+1)}\Biggr \} \cdots {\tilde {U}}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}\Biggr \} {\tilde {U}}_{{t}_{l},0}^{(k+l)}. \end{eqnarray} \tag{ A.3 }$$

Then B₁ is expressed as

$$\begin{eqnarray} &&{B}_{1}=\sum _{l=0}^{L-1}\left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k+l}{B}_{1}^{(k+l)}{Q}_{N,0}^{(k+l)}\vphantom {\sum _{1}^{2}} \right \vert . \end{eqnarray} \tag{ A.4 }$$

In equation (A.3), a quantity like

$$\begin{eqnarray*} &&\frac{{\partial }^{n}{\tilde {U}}_{t,s}^{(k)}}{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})\partial {\alpha }_{{\sigma }_{2}}({t}_{2})\cdots \partial {\alpha }_{{\sigma }_{n}}({t}_{n})} \end{eqnarray*}$$

appears. We find that

$$\begin{eqnarray} \frac{\partial {\tilde {U}}_{t,s}^{(k)}}{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}&=&-\mathrm{i}\int \nolimits _{s}^{t}\mathrm{d}{t}^{\prime}\hspace{0.167em} {\tilde {U}}_{t,{t}^{\prime}}^{(k)}\frac{\partial {\mathcal{L}}_{a}(\alpha ({t}^{\prime}))}{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}{\tilde {U}}_{{t}^{\prime},s}^{(k)}\nonumber\\ &=&-\mathrm{i}g\int \nolimits _{s}^{t}\mathrm{d}{t}^{\prime}\hspace{0.167em} {\tilde {U}}_{t,{t}^{\prime}}^{(k)}\frac{\partial {\alpha }_{{\sigma }_{1}}({t}^{\prime})}{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}\sum _{j=1}^{k}{X}_{j}^{(-{\sigma }_{1})\times }{\tilde {U}}_{{t}^{\prime},s}^{(k)}\nonumber\\ &=&-\mathrm{i}g\int \nolimits _{s}^{t}\mathrm{d}{t}^{\prime}\hspace{0.167em} {\tilde {U}}_{t,{t}^{\prime}}^{(k)}{g}_{a a}({t}^{\prime}-{t}_{1})\sum _{j=1}^{k}{X}_{j}^{(-{\sigma }_{1})\times }{\tilde {U}}_{{t}^{\prime},s}^{(k)}, \end{eqnarray} \tag{ A.5 }$$

where ${X}_{j}^{(-{\sigma }_{1})\times }(\cdot ):= [{X}_{j}^{(-{\sigma }_{1})},(\cdot )]$ is the commutator. Remember that |g_aa(t)| ≤ 1 because of equation (48). Therefore,

$$\begin{eqnarray} &&\left \Vert \frac{\partial {\tilde {U}}_{t,s}^{(k)}}{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}\right \Vert _{\infty }\leq 2 k g\Vert X\Vert _{\infty }(t-s). \end{eqnarray} \tag{ A.6 }$$

Similarly, we can obtain

$$\begin{eqnarray} &&\left \Vert \frac{{\partial }^{n}{\tilde {U}}_{t,s}^{(k)}}{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})\partial {\alpha }_{{\sigma }_{2}}({t}_{2})\cdots \partial {\alpha }_{{\sigma }_{n}}({t}_{n})}\right \Vert _{\infty }\leq \left [2 k g\Vert X\Vert _{\infty }(t-s)\right ]^{n}. \end{eqnarray} \tag{ A.7 }$$

Because of equation (A.7), we can replace $\partial {\tilde {U}}_{t,s}/\partial {\alpha }_{\sigma }({t}^{\prime})$ simply by $2 k g\Vert X\Vert _{\infty }(t-s){\tilde {U}}_{t,s}$ as far as it is concerned with the upper bound of the operator norm. This fact makes it much easier to evaluate the upper bound.

We define the set

$$\begin{eqnarray*} &&{\mathcal{N}}_{l,n}:= \left \{ \{ {m}_{1},{m}_{2},\ldots ,{m}_{n}\} :1\leq {m}_{1}\lt {m}_{2}\lt \cdots \lt {m}_{n}\leq l,\quad {m}_{j}\in \mathbb{N}\right \} . \end{eqnarray*}$$

If one of the elements of _l,n is denoted by N_l,n = {m₁,m₂,...,m_n}, we define ${N}_{l,n}^{c}:= \{ 1,2,\ldots ,l\} \setminus {N}_{l,n}$.

By using equation (A.7), we then obtain

$$\begin{eqnarray} &&\fl \left \Vert {B}_{1}^{(k+l)}\right \Vert _{\infty }\leq \frac{(g\Vert X\Vert _{\infty }t)^{l}}{l !}\mathop{\sum }\limits _{\{ {\sigma }_{j}=\pm 1\} }\sum _{n=0}^{l}\mathop{\sum }\limits _{{N}_{l,n}\in {\mathcal{N}}_{l,n}}\nonumber\\ &&\times ~\left \Vert \left (\mathop{\prod }\limits _{j\in {N}_{l,n}}\frac{\partial }{\partial {\alpha }_{{\sigma }_{j}}(t)}\right ){\mathcal{O}}^{(k)}(\alpha (t),{\alpha }^{\ast }(t),\{ {\beta }_{r}(t),\beta _{r}^{\ast }(t)\} )\right \Vert _{\infty }\nonumber\\ &&\times ~(k+l-1)(k+l-2)\cdots (k+n)\times (2 g\Vert X\Vert _{\infty }t)^{l-n}. \end{eqnarray} \tag{ A.8 }$$

Because ^(k)∈ℬ, we have

$$\begin{eqnarray} &&\fl \left \Vert {B}_{1}^{(k+l)}\right \Vert _{\infty }\leq \frac{(g\Vert X\Vert _{\infty }t)^{l}}{l !}\mathop{\sum }\limits _{\{ {\sigma }_{j}=\pm 1\} }\sum _{n=0}^{l}\mathop{\sum }\limits _{{N}_{l,n}\in {\mathcal{N}}_{l,n}}\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{\kappa }^{n}\nonumber\\ &&\times ~(k+l-1)(k+l-2)\cdots (k+n)\times (2 g\Vert X\Vert _{\infty }t)^{l-n}. \end{eqnarray} \tag{ A.9 }$$

Since ∑_{σj=±1}1 = 2^l, ∑_Nl,n∈l,n1 = l!/(n!(l − n)!), we obtain the upper bound

$$\begin{eqnarray} &&\fl \left \Vert {B}_{1}^{(k+l)}\right \Vert _{\infty }\leq \sum _{n=0}^{\infty }\frac{{\kappa }^{n}}{n !}\frac{(k+l-1)!}{(l-n)!(k+n-1)!}(2 g\Vert X\Vert _{\infty }t)^{2 l-n}\nonumber\\ &&\times ~\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}. \end{eqnarray} \tag{ A.10 }$$

Here, we used (k + l − 1)!/[(l − n)!(k + n − 1)!] ≤ 2^k+l−1 and restricted the time interval t∈[0,t₀] where t₀ is determined by 2g||X||_∞t₀ = 1/2. This restriction will be removed in section 5.2. Then we obtain

$$\begin{eqnarray} &&\left \Vert {B}_{1}^{(k+l)}\right \Vert _{\infty }\leq {2}^{k-l-1}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}. \end{eqnarray} \tag{ A.11 }$$

From equations (A.4) and (A.11), we obtain

$$\begin{eqnarray} &&\fl {B}_{1}\leq {2}^{k}{\mathrm{e}}^{2 \kappa }\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,0}^{(0)}. \end{eqnarray} \tag{ A.12 }$$

Here, |c₀+c_aα(t)+∑_rc_rβ_r(t)|^q∈ℬ (which will be explicitly shown in section 5.1). By assumption, at the initial time t = 0, the expectation value is finite for any ^(k)∈ℬ. Therefore, B₁ is also finite.

A.2. Evaluation of B₂

B₂ is given by

$$\begin{eqnarray} &&\fl {B}_{2}=\Biggl \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\nonumber\\ &&\tqs \times ~{U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {U}_{{t}_{L-1},{t}_{L}}^{(k+L-1)}{\mathcal{W}}^{(k+L-1)}{Q}_{N,{t}_{L}}^{(k+L)}\Biggr \vert . \end{eqnarray} \tag{ A.13 }$$

This is written in the form

$$\begin{eqnarray} &&\fl {B}_{2}\leq \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\hspace{0.167em} {\mathrm{Tr}}_{1,2,\ldots ,k+L}\left \Vert {B}_{2}^{(k+L)}\right \Vert _{\infty }{Q}_{N,{t}_{L}}^{(k+L)},\nonumber\\ \end{eqnarray} \tag{ A.14 }$$

where ${B}_{2}^{(k+L)}$ is given by

$$\begin{eqnarray} &&\fl {B}_{2}^{(k+L)}:= {g}^{L}\sum _{\{ {\sigma }_{j}=\pm 1\} }^{l}{X}_{k+1}^{({\sigma }_{1})}{X}_{k+2}^{({\sigma }_{2})}\cdots {X}_{k+L}^{({\sigma }_{L})}\frac{\partial }{\partial {\alpha }_{{\sigma }_{L}}({t}_{L})}\Biggl \{ \frac{\partial }{\partial {\alpha }_{{\sigma }_{L-1}}({t}_{L-1})}\Biggl \{ \cdots \frac{\partial }{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}\nonumber\\ &&\times ~\Biggl \{ {\mathcal{O}}^{(k)}\left (\alpha (t),{\alpha }^{\ast }(t),\{ {\beta }_{r}(t),{\beta }_{r}^{\ast }(t)\} \right ){\tilde {U}}_{t,{t}_{1}}^{(k)}\Biggr \} {\tilde {U}}_{{t}_{1},{t}_{2}}^{(k+1)}\Biggr \} \cdots {\tilde {U}}_{{t}_{L-1},{t}_{L}}^{(k+L-1)}\Biggr \} . \end{eqnarray} \tag{ A.15 }$$

Similarly to the evaluation of B₁, we obtain the upper bound

$$\begin{eqnarray} &&\fl {B}_{2}\leq \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\cdots \int \nolimits _{0}^{{t}_{L}-1}\mathrm{d}{t}_{L}\hspace{0.167em} L !(2 g\Vert X\Vert _{\infty })^{L}\nonumber\\ &&\times \sum _{m=0}^{L}\frac{1}{m !}{2}^{k+L-1}(2 g\Vert X\Vert _{\infty }t)^{L-m}{\kappa }^{m}\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t) + \mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}_{L}}^{(k+L)}\nonumber\\ &&\leq ~(2 g\Vert X\Vert _{\infty }t)^{L}\sum _{m=0}^{L}\frac{1}{m !}{2}^{k+L-1}(2 g\Vert X\Vert _{\infty }t)^{L-m}\nonumber\\ &&\times ~{\kappa }^{m}\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}^{\prime}}^{(k+L)}.\nonumber\\ \end{eqnarray} \tag{ A.16 }$$

Here, we again restrict the time interval t∈[0,t₀] (t₀ was determined by 2g||X||_∞t₀ = 1/2 in section A.1). Then we obtain

$$\begin{eqnarray} &&\fl \left \Vert {B}_{2}^{(k+L)}\right \Vert _{\infty }\leq {2}^{-L+k-1}\sum _{m=0}^{L}\frac{1}{m !}(2 \kappa )^{m}\mathop{\sum }\limits _{q=0}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}\nonumber\\ &&\leq ~{2}^{-L+k-1}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}. \end{eqnarray} \tag{ A.17 }$$

Substituting it into equation (A.14), we have

$$\begin{eqnarray} &&\fl {B}_{2}\leq {2}^{-L+k-1}\sum _{m=0}^{L}\frac{1}{m !}(2 \kappa )^{m}\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\nonumber\\ &&\times ~\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}^{\prime}}^{(0)}\nonumber\\ &&\leq ~{2}^{-L+k-1}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\nonumber\\ &&\times ~\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}^{\prime}}^{(0)}. \end{eqnarray} \tag{ A.18 }$$

A.3. Evaluation of B₃

Next we evaluate B₃, which is given by

$$\begin{eqnarray} &&\fl {B}_{3}=\frac{1}{N}\sum _{l=1}^{L}\left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\right .\nonumber\\ &&\times ~\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\hspace{0.167em} {U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\nonumber\\ &&\left .\cdots {U}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}{\mathcal{V}}^{(k+l-1)}{Q}_{N,{t}_{l}}^{(k+l-1)}\vphantom {\sum _{1}^{2}} \right \vert . \end{eqnarray} \tag{ A.19 }$$

If we put

$$\begin{eqnarray} &&\fl {B}_{3}^{(k+l)}:= (-\mathrm{i}g)^{l}t\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\nonumber\\ &&\times ~\mathop{\sum }\limits _{\{ {\sigma }_{j}=\pm 1\} }X_{k+1}^{({\sigma }_{1})}{X}_{k+2}^{({\sigma }_{2})}\cdots {X}_{k+l}^{({\sigma }_{l})}\nonumber\\ &&\times ~\frac{\partial }{\partial {\alpha }_{{\sigma }_{l}}({t}_{l})}\biggl \{ \frac{\partial }{\partial {\alpha }_{{\sigma }_{l-1}}({t}_{l-1})}\biggl \{ \cdots \frac{\partial }{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}\biggl \{ \nonumber\\ &&\times ~{\mathcal{O}}^{(k)}(\alpha (t),{\alpha }^{\ast }(t),\{ {\beta }_{r}(t),{\beta }_{r}^{\ast }(t)\} ){\tilde {U}}_{t,{t}_{1}}^{(k)}\biggr \} {\tilde {U}}_{{t}_{1},{t}_{2}}^{(k+1)}\biggr \} \cdots {\tilde {U}}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}\biggr \} {\mathcal{V}}_{{\sigma }_{l}}^{(k+l)}, \end{eqnarray} \tag{ A.20 }$$

where

$$\begin{eqnarray} {\mathcal{V}}_{\sigma }^{(k)}(\cdot ):= \left \{ \begin{array}{@{}ll@{}} \displaystyle (\cdot )\sum _{j=1}^{k}{X}_{j}\quad &\displaystyle \text{for }\sigma =1,\\ \displaystyle \sum _{j=1}^{k}{X}_{j}^{\dagger }(\cdot )\quad &\displaystyle \text{for }\sigma =-1, \end{array}\right . \end{eqnarray} \tag{ A.21 }$$

we obtain the upper bound of B₃,

$$\begin{eqnarray} &&{B}_{3}\leq \frac{1}{N}\sum _{l=1}^{L}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\left \Vert {B}_{3}^{(k+l)}\right \Vert _{\infty }{Q}_{N,{t}^{\prime}}^{(0)}. \end{eqnarray} \tag{ A.22 }$$

Similarly to the evaluation of B₂, we obtain for 0 < t ≤ t₀

$$\begin{eqnarray} &&\left \Vert {B}_{3}^{(k+l)}\right \Vert _{\infty }\leq (k+l){2}^{k-l-1}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}. \end{eqnarray} \tag{ A.23 }$$

Therefore,

$$\begin{eqnarray} &&\fl {B}_{3}\leq \frac{1}{N}\sum _{l=1}^{L}(k+l){2}^{k-l-1}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\nonumber\\ &&\times ~\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}_{l}}^{(0)}. \end{eqnarray} \tag{ A.24 }$$

By elementary calculation, we find ${\mathop{\sum }\nolimits }_{l=1}^{L}l/{2}^{l}\lt 2$, and thus

$$\begin{eqnarray} &&\fl {B}_{3}\lt \frac{1}{N}{2}^{k-1}(k+2){\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}^{\prime}}^{(0)}.\nonumber\\ \end{eqnarray} \tag{ A.25 }$$

A.4. Evaluation of B₄

Finally, we evaluate the contribution of B₄. Its upper bound is given by

$$\begin{eqnarray} &&{B}_{4}\leq \frac{1}{N}\sum _{l=1}^{L}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\left \Vert {B}_{4}^{(k+l)}\right \Vert _{\infty }{Q}_{N,{t}^{\prime}}^{(k+l)}, \end{eqnarray} \tag{ A.26 }$$

where ${B}_{4}^{(k+l)}$ is defined by

$$\begin{eqnarray} &&\fl {B}_{4}^{(k+l)}:= (-\mathrm{i}g)^{l}(k+l-1)\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\nonumber\\ &&\times ~\mathop{\sum }\limits _{\{ {\sigma }_{j}=\pm 1\} }X_{k+1}^{({\sigma }_{1})}{X}_{k+2}^{({\sigma }_{2})}\cdots {X}_{k+l}^{({\sigma }_{l})}\nonumber\\ &&\times ~\frac{\partial }{\partial {\alpha }_{{\sigma }_{l}}({t}_{l})}\left \{ \frac{\partial }{\partial {\alpha }_{{\sigma }_{l-1}}({t}_{l-1})}\left \{ \cdots \frac{\partial }{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}\Biggl \{ \right .\right .\nonumber\\ &&\left .\left .\vphantom {\frac {\partial }{\partial {\alpha }_{{\sigma }_{l}}({t}_{l})}\{ \frac {\partial }{\partial {\alpha }_{{\sigma }_{l-1}}({t}_{l-1})} \cdots \frac {\partial }{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}} \times ~{\mathcal{O}}^{(k)}(\alpha (t),{\alpha }^{\ast }(t),\{ {\beta }_{r}(t),{\beta }_{r}^{\ast }(t)\} ){\tilde {U}}_{t,{t}_{1}}^{(k)}\Biggr \} {\tilde {U}}_{{t}_{1},{t}_{2}}^{(k+1)}\right \} \cdots {\tilde {U}}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}\right \} . \end{eqnarray} \tag{ A.27 }$$

We can evaluate the upper bound of ${B}_{4}^{(k+l)}$ by some calculations similar to those in previous subsections, and the result is

$$\begin{eqnarray} &&\left \Vert {B}_{4}^{(k+l)}\right \Vert _{\infty }\leq (k+l-1){2}^{k-l-1}{\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}, \end{eqnarray} \tag{ A.28 }$$

for 0 ≤ t ≤ t₀. Therefore, we conclude that

$$\begin{eqnarray} &&\fl {B}_{4}\lt \frac{1}{N}{2}^{k-1}(k+1){\mathrm{e}}^{2 \kappa }\sum _{q=0}^{s}{a}_{q}\mathop{\max }\limits _{{t}^{\prime}\in [0,t]}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right )\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}{Q}_{N,{t}^{\prime}}^{(0)},\nonumber\\ \end{eqnarray} \tag{ A.29 }$$

in the time interval 0 ≤ t ≤ t₀.

B.1. Evaluation of ${B}_{1}^{\prime}$

First, we consider the term ${B}_{1}^{\prime}$. Similarly to the case of B₁, we obtain the following expression:

$$\begin{eqnarray} &&\fl {B}_{1}^{\prime}=\sum _{l=0}^{L-1}\Biggl \vert \int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{l-1}}\mathrm{d}{t}_{l}\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k+l}\nonumber\\ &&\times ~\sum _{\{ {\sigma }_{j}=\pm 1\} }^{l}{X}_{k+1}^{({\sigma }_{1})}{X}_{k+2}^{({\sigma }_{2})}\cdots {X}_{k+l}^{({\sigma }_{l})}\nonumber\\ &&\times ~\frac{\partial }{\partial {\alpha }_{{\sigma }_{l}}({t}_{l})}\biggl \{ \frac{\partial }{\partial {\alpha }_{{\sigma }_{l-1}}({t}_{l-1})}\biggl \{ \cdots \frac{\partial }{\partial {\alpha }_{{\sigma }_{1}}({t}_{1})}\biggl \{ \nonumber\\ &&\times ~{\mathcal{O}}^{(k)}(\alpha (t),{\alpha }^{\ast }(t)){\tilde {U}}_{t,{t}_{1}}^{(k)}\biggr \} {\tilde {U}}_{{t}_{1},{t}_{2}}^{(k+1)}\biggr \} \cdots {\tilde {U}}_{{t}_{l-1},{t}_{l}}^{(k+l-1)}\biggr \} {\tilde {U}}_{{t}_{l},0}^{(k+l)}\left ({Q}_{N,0}^{(k+l)}-{Q}_{\mathrm{ MF},0}^{(k+l)}\right )\Biggr \vert .\nonumber\\ \end{eqnarray} \tag{ B.1 }$$

It is expressed as

$$\begin{eqnarray} &&{B}_{1}^{\prime}=\sum _{l=0}^{L-1}\left \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k+l}{B}_{1}^{(k+l)}\left ({Q}_{N,0}^{(k+l)}-{Q}_{\mathrm{ MF},0}^{(k+l)}\right )\right \vert , \end{eqnarray} \tag{ B.2 }$$

where ${B}_{1}^{(k+l)}$ was defined by equation (A.3). From the assumption of the initial condition (49), we can conclude that ${\lim }_{N\rightarrow \infty }{B}_{1}^{\prime}=0$ if we can show ${B}_{1}^{(k+l)}\in \mathcal{B}$.

In order to show ${B}_{1}^{(k+l)}\in \mathcal{B}$, we must consider the quantity

$$\begin{eqnarray} &&\left \Vert \frac{{\partial }^{m}{B}_{1}^{(k+l)}}{\partial {\alpha }_{{\tau }_{1}}\partial {\alpha }_{{\tau }_{2}}\cdots \partial {\alpha }_{{\tau }_{m}}}\right \Vert _{\infty }, \end{eqnarray} \tag{ B.3 }$$

where τ_j =± 1.

A calculation analogous to that for B₁ yields

$$\begin{eqnarray} &&\fl \left \Vert \frac{{\partial }^{m}{B}_{1}^{(k+l)}}{\partial {\alpha }_{{\tau }_{1}}\partial {\alpha }_{{\tau }_{2}}\cdots \partial {\alpha }_{{\tau }_{m}}}\right \Vert _{\infty }\leq \frac{(g\Vert X\Vert _{\infty }t)^{l}}{l !}\mathop{\sum }\limits _{\{ {\sigma }_{j}=\pm 1\} }\sum _{n=0}^{l}\mathop{\sum }\limits _{{N}_{l,n}\in {\mathcal{N}}_{l,n}}\sum _{p=0}^{m}\mathop{\sum }\limits _{{N}_{m,p}\in {\mathcal{N}}_{m,p}}\nonumber\\ &&\times ~\left \Vert \left (\mathop{\prod }\limits _{j\in {N}_{l,n}}\frac{\partial }{\partial {\alpha }_{{\sigma }_{j}}(t)}\right )\left (\mathop{\prod }\limits _{j\in N_{m,p}^{c}}\frac{\partial }{\partial {\alpha }_{{\tau }_{j}}(t)}\right ){\mathcal{O}}^{(k)}(\alpha (t),{\alpha }^{\ast }(t))\right \Vert _{\infty }\nonumber\\ &&\times ~(k+l)^{p}(k+l-1)(k+l-2)\cdots (k+n)\times (2 g\Vert X\Vert _{\infty }t)^{l-n+p}. \end{eqnarray} \tag{ B.4 }$$

Because ^(k)∈ℬ, we have

$$\begin{eqnarray} &&\fl \left \Vert \frac{{\partial }^{m}{B}_{1}^{(k+l)}}{\partial {\alpha }_{{\tau }_{1}}\partial {\alpha }_{{\tau }_{2}}\cdots \partial {\alpha }_{{\tau }_{m}}}\right \Vert _{\infty }\leq \frac{(g\Vert X\Vert _{\infty }t)^{l}}{l !}\mathop{\sum }\limits _{\{ {\sigma }_{j}=\pm 1\} }\sum _{n=0}^{l}\mathop{\sum }\limits _{{N}_{l,n}\in {\mathcal{N}}_{l,n}}\sum _{p=0}^{m}\mathop{\sum }\limits _{{N}_{m,p}\in {\mathcal{N}}_{m,p}}\nonumber\\ &&\times ~(k+l)^{p}(k+l-1)(k+l-2)\cdots (k+n)(2 g\Vert X\Vert _{\infty }t)^{l-n+p}\nonumber\\ &&\times ~{\kappa }^{n-m+p}\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}\right \vert ^{q}. \end{eqnarray} \tag{ B.5 }$$

Since ∑_{σj=±1}1 = 2^l, ∑_Nl,n∈l,n1 = l!/(n!(l − n)!), we obtain the upper bound

$$\begin{eqnarray} &&\fl \left \Vert \frac{{\partial }^{m}{B}_{1}^{(k+l)}}{\partial {\alpha }_{{\tau }_{1}}\partial {\alpha }_{{\tau }_{2}}\cdots \partial {\alpha }_{{\tau }_{m}}}\right \Vert _{\infty }\leq {2}^{k+l-1}\sum _{n=0}^{l}\frac{1}{n !}{\kappa }^{n}(2 g\Vert X\Vert _{\infty }t)^{2 n-l}\nonumber\\ &&\times ~\left [\kappa +2(k+l)g\Vert X\Vert _{\infty }t\right ]^{m}\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}. \end{eqnarray} \tag{ B.6 }$$

From this expression, we find that ${B}_{1}^{(k+l)}\in \mathcal{B}$ for any fixed t. Thus we can conclude that

$$\begin{eqnarray} &&\mathop{\lim }\limits _{N\rightarrow \infty }B_{1}^{\prime}=0. \end{eqnarray} \tag{ B.7 }$$

B.2. Evaluation of ${B}_{2}^{\prime}$

${B}_{2}^{\prime}$ is given by

$$\begin{eqnarray} &&\fl {B}_{2}^{\prime}=\Biggl \vert \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}{\mathcal{O}}^{(k)}\left (\alpha ,{\alpha }^{\ast },\{ {\beta }_{r},{\beta }_{r}^{\ast }\} \right )\int \nolimits _{0}^{t}\mathrm{d}{t}_{1}\int \nolimits _{0}^{{t}_{1}}\mathrm{d}{t}_{2}\cdots \int \nolimits _{0}^{{t}_{L-1}}\mathrm{d}{t}_{L}\nonumber\\ &&\times ~{U}_{t,{t}_{1}}^{(k)}{\mathcal{W}}^{(k)}{U}_{{t}_{1},{t}_{2}}^{(k+1)}{\mathcal{W}}^{(k+1)}\cdots {U}_{{t}_{L-1},{t}_{L}}^{(k+L-1)}{\mathcal{W}}^{(k+L-1)}\left ({Q}_{N,{t}_{L}}^{(k+L)}-{Q}_{\mathrm{ MF},{t}_{L}}^{(k+L)}\right )\Biggr \vert . \end{eqnarray} \tag{ B.8 }$$

This is written in the form

$$\begin{eqnarray} &&{B}_{2}^{\prime}\leq \int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}\left \Vert {B}_{2}^{(k+L)}\right \Vert _{\infty }\left ({Q}_{N,{t}_{L}}^{(k+L)}+{Q}_{\mathrm{ MF},{t}_{L}}^{(k+L)}\right ), \end{eqnarray} \tag{ B.9 }$$

where ${B}_{2}^{(k+L)}$ is given by equation (A.15). From the upper bound of ${B}_{2}^{(k+L)}$ (see equation (A.17)), it yields

$$\begin{eqnarray} &&\fl {B}_{2}^{\prime}\leq {2}^{-L+k-1}{\mathrm{e}}^{2 \kappa }\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}\nonumber\\ &&\times \left ({Q}_{N,{t}_{L}}^{(k+L)}+{Q}_{\mathrm{ MF},{t}_{L}}^{(k+L)}\right ), \end{eqnarray} \tag{ B.10 }$$

which is almost the same as equation (A.18). Because

$$\begin{eqnarray*} &&\fl \mathop{\lim }\limits _{N\rightarrow \infty }\int {\mathrm{d}}^{2}\alpha \left (\mathop{\prod }\limits _{r}\int {\mathrm{d}}^{2}{\beta }_{r}\right ){\mathrm{Tr}}_{1,2,\ldots ,k}\sum _{q=0}^{s}{a}_{q}\left \vert {c}_{0}+{c}_{a}\alpha (t)+\mathop{\sum }\limits _{r}{c}_{r}{\beta }_{r}(t)\right \vert ^{q}\left ({Q}_{N,{t}_{L}}^{(k+L)}+{Q}_{\mathrm{ MF},{t}_{L}}^{(k+L)}\right ) \end{eqnarray*}$$

is finite and not diverging as L → ∞, it is concluded that

$$\begin{eqnarray} &&\mathop{\lim }\limits _{L\rightarrow \infty }\mathop{\lim }\limits _{N\rightarrow \infty }B_{2}^{\prime}=0, \end{eqnarray} \tag{ B.11 }$$

for 0 ≤ t ≤ t₀.