Quantum input-output theory : Why do we multiply by density of mode to have a number of photon **per unit of time** In this paper, https://journals.aps.org/pra/abstract/10.1103/PhysRevA.31.3761, we work with input-output theory. I will first summarize the physics of it and then ask my question.
In input-output theory we model everything by saying : I send an input field, it interacts with a quantum system, and after the interaction I have an output field.
==Input==> [Interaction of field with Q.Syst] ==Output==> 


*

*The Hamiltonian is :
$$ H = H_s + H_{field} + H_{int} $$
Where $H_{field}=\int d\omega  ~ \hbar \omega b^{\dagger}(\omega)
   b(\omega) $. It is not necessary to explicit the other parts for my question

*In Heisenberg picture, if we didn't had any interaction, the field
would freely evolve such that 
$$b(\omega,t)=e^{-i \omega t} b^0(\omega)$$

Before moving further, it is important to underline that $b^{\dagger}(\omega)
   b(\omega)$ is not the number of photon of frequency $\omega$ because it has the dimension of a time. 
This is basically what will cause me all the problem. So for me if I work with a discrete number of mode, I would write :
$$ \sum_k \hbar \omega_k a^{\dagger}(\omega_k) a(\omega_k) = \int d \omega ~ \hbar \omega \nu(\omega) a^{\dagger}(\omega) a(\omega) $$ and I would identify $b(\omega) = \sqrt{\nu (\omega)} a(\omega)$ where $\nu(\omega) d \omega$ is the number of mode I have in $[\omega; \omega+ d \omega]$.



*

*We define the input field as :
$$b_{in}(t)=\frac{1}{\sqrt{2 \pi}} \int d\omega ~ e^{-i \omega t} 
   b^0(\omega)$$
It looks like a Fourier transform but the way I understand it is more : we make evolve all mode at time $t$ in Heisenberg picture assuming they are not interacting (which is the case of the input field before the interaction), and we sum on those modes : it is the definition of the total input field.
My question :
In this paper (and more generally everytime input-output theory is used), they say that $<b^{\dagger}_{in}(t) b_{in}(t)>$ is the number of photon per unit of time I have at time $t$ in my input field.
I don't understand this. I agree it has the good dimension but why would this quantity physically represent it ?
 A: Some remarks:


*

*From the question (emphasis mine)

Before moving further, it is important to underline that $b^†(ω)b(ω)$ is not the number of photon of frequency $ω$ because it has the dimension of the inverse of a time.

How so? According to 

$H_{field}=\int d\omega  ~ \hbar \omega b^{\dagger}(\omega)
   b(\omega)$

we ought to have that $ \int d\omega  ~ \hbar \omega \langle b^{\dagger}(\omega)
   b(\omega)\rangle$ is the total energy of the bath in the absence of interactions with the system. Therefore $\langle b^{\dagger}(\omega)
   b(\omega)\rangle$ is simply the number of photons per unit frequency. The dimension is therefore inverse frequency, not inverse time.

*
$\sum_k \hbar \omega_k a^{\dagger}(\omega_k) a(\omega_k) = \int d \omega ~ \hbar \omega \nu(\omega) a^{\dagger}(\omega) a(\omega)$

It is worth noting here that in this step an approximation is being performed, namely the discrete system modes are being replaced by a continuum. Also note that the notation is slightly abusive, since the $a^\dagger$ operator changes units from the left to the right.

*
We define the input field as :
$b_{in}(t)=\frac{1}{\sqrt{2 \pi}} \int d\omega ~ e^{-i \omega t} 
   b^0(\omega)$
It looks like a Fourier transform but the way I understand it is more : we make evolve all mode at time t
  in Heisenberg picture assuming they are not interacting (which is the case of the input field before the interaction), and we sum on those modes : it is the definition of the total input field.

I quite like this interpretation. Let me rephrase it a bit: at the start you have a quantum mechanical "wave packet" that is built from a superposition of bath modes, each evolving freely at frequency $\omega$. Note that this contains the notion of being an boundary condition at time $t_0$ (indicated only by the zero superscript in this formula). In many cases this will be considered asymptotically, with $t_0$ approaching the infinite past.

*
$b_{in}(t)=\frac{1}{\sqrt{2 \pi}} e^{-i \omega t} 
   b^0(\omega)=\frac{1}{\sqrt{2 \pi}} e^{-i \omega t} \sqrt{\nu(\omega)} a^0(\omega)$

This formula in the question is not completely correct, especially the continuum system modes on the right hand side. It is certainly not what Gardiner&Collett have in their paper (see formula (2.22)). They only have a single discrete mode. If you want to have a continuum of system modes, there should at least be an integral for that somewhere, unless your coupling is local in frequency. But the latter would just correspond to a single mode problem with messy notation again. Either way the units are wrong in this, as pointed out in remark 2.
I initially thought that this was where the confusion came from, but after StarBucK's I am adding this edit to address the real question:
EDIT:
So now that we have understood the units of $b(\omega)$, which was also nicely expained again in an answer by jgerber that was posted since, we can look at the units of $b_\textrm{in}(t)$. To understand this let us investigate the definition of the input operators a bit further. The original bath operator $b(\omega)$ is in the Heisenberg picture (as also pointed out by jgerber). So this operator is already time dependent and could (or maybe should) be written $b(t, \omega)$. As we saw above, what $\langle b^\dagger(t, \omega) b(t, \omega) \rangle$ then represents physically is the number of photons per unit frequency (so "per mode") at time t (not per unit time). So in other words: $b(t, \omega)$ is our standard photon operator, just for a continuum, not for a discrete mode.
The definition of the input operator can then be written:
$$b_{\textrm{in}}(t)=\frac{1}{\sqrt{2 \pi}} \int d\omega ~ e^{-i \omega (t-t_0)} b(t_0, \omega)$$
Note that $e^{i \omega t_0} b(t_0, \omega)$ is physically the photon operator in the interaction picture (that is with the free time evolution taken out) at time $t_0$. So if you have no interactions, then $e^{i \omega t_0} b(t_0, \omega)$ is actually independent of $t_0$. This means $b_{\textrm{in}}(t)$ is really the Fourier transform of the interaction picture operator at time $t_0$.
To make it a bit clearer what I am saying, you can also define an input operator in the frequency domain. The definition is just a Fourier transform again and if we evaluate this Fourier transform, we get a very simple result:
$$b_{\textrm{in}}(\omega) = \frac{1}{\sqrt{2 \pi}} \int dt e^{i \omega t} b_{\textrm{in}}(t) = e^{i \omega t_0} b(t_0, \omega).$$
So the input operator in the frequency domain is exactly the interaction picture photon operator at time $t_0$! Mathematically this is all simple, we are just doing Fourier transforms back and forth. But physically, this gives a lot of insight into what the input operators mean in my opinion.
So say we say that the expectation value of these frequency space input operators is some function $\langle b^\dagger_{\textrm{in}}(\omega) b_{\textrm{in}}(\omega) \rangle = I(\omega)$. I have called this function $I$ on purpose, because this represents the intensity spectrum that you send into your system at time $t_0$. So if you have some wavepacket flying towards your cavity/interaction region, the input operators give you the spectrum of this wavepacket.
The time-frequency relation then behaves very similarly to classical optics. We have
$$ \omega\textrm{-domain amplitude} \xleftarrow[]{\textrm{Expectation value}}b_{\textrm{in}}(\omega) \xrightarrow[]{\textrm{Fourier transform}}  b_{\textrm{in}}(t) \xrightarrow[]{\textrm{Expectation value}} t\textrm{-domain amplitude}$$
So my advice is to think about it in terms of wavepackets. $\langle b^\dagger_{\textrm{in}}(\omega) b_{\textrm{in}}(\omega) \rangle$ gives you the number of photons per unit frequency at frequency $\omega$ in the wavepacket. $\langle b^\dagger_{\textrm{in}}(t) b_{\textrm{in}}(t) \rangle$ gives you the number of photons per unit time at time $t$ in the wavepacket. Here is a picture (source): 

The picture is for classical fields, so you do not have the whole business with expectation values, but the principle is the same. $|E(\omega)|^2$ is the intensity per unit frequency at frequency $\omega$, $|E(t)|^2$ is the intensity per unit time at time t.
Summary: after stripping away the weirdness of the definition of input operators and the interaction picture, this is really just Fourier transforming wave packets.
A: Disclaimer: I struggle with this question myself a bit. I think I have something to add that may be helpful to you but it isn't going to be a superpolished answer and it might be a little unclear at times.
Consider
$$
H = \int \hbar \omega b^{\dagger}(\omega)b(\omega) d\omega
$$
We see that, since $H$ has units of energy, $\hbar \omega$ has units of energy, and $d\omega$ has units of frequency we can see that $b^{\dagger}(\omega)b(\omega)$ has units of inverse frequency, that is $Hz^{-1}$.
In this form it is clear that $\hbar \omega b^{\dagger}(\omega)b(\omega)$ represents the total energy per unit frequency. That is, we can see that $\hbar \omega b^{\dagger}(\omega)b(\omega)$ is like a power spectral density for the photon field. If you integrate over some frequency range $[\omega_1, \omega_2]$ you get the total energy contained in that bandwidth. $b^{\dagger}(\omega)b(\omega)$ is then some sort of excitation spectral density. Again, if you integrate you get the total number of excitations in some bandwidth. If $b^{\dagger}(\omega)b(\omega)$ is the power spectral density of the signal then $b(\omega)$ is the signal itself, it is related to the square root of the power spectral density. Thus it has units of $Hz^{-\frac{1}{2}}$. This unit comes up often in signal processing. We talk about "Something per root Hz" when we want to discuss a power spectral density which is scaled to have units comparable to the amplitude of a signal rather than the power in the signal.
Ok, so to recap $b(\omega)$ has units of $Hz^{-\frac{1}{2}}$. We now consider the "in" field:
$$
b_{in}(t) = \int e^{i\omega t}b_0(\omega) d\omega
$$
First notice that $b_{in}(t)$ is a function of time not frequency. This is because frequency was integrated over. Next, you've noticed that $b_{in}(t)$ has units of $Hz^{\frac{1}{2}}$ now, or $s^{-\frac{1}{2}}$. So I think this should first and foremost simply be thought of as an integral over bath operators. In the description with $b_0(\omega)$ we have split thing up by their frequency modes, somehow we've already done a Fourier transform. $b_{in}(t)$ captures the integrated signal which is the composition of all of those modes into one operator.
But there is an important point which has not been discussed yet. The input-output formalism is often considering photons inside of a cavity described by some cavity operators $a$ which are coupled to some continuum modes outside of the cavity, $b(\omega)$. Let's think about the electromagnetic fields inside and outside of the cavity. 
Inside the cavity there are standing wave modes of light. Note that the standing wave is composed of a left travelling wave and right travelling wave. There is no spatial translation of the energy in the cavity. You can think of the cavity as a box and the total energy in the box is $\hbar \omega a^{\dagger} a$ (no funny integral, now $a$ is dimensionless in contrast to $b(\omega)$ which has units of $Hz^{-\frac{1}{2}}$ or $b_{in}(t)$ which has units of $s^{-\frac{1}{2}}$).
However, in free space it is possible to have a travelling wave in one direction without also having the travelling wave in the other direction. That is, it is possible to have waves carrying energy in a particular direction, we can talk now about energy flow.
Both in a cavity and in free space the electromagnetic field must satisfy the Helmholtz equation. The difference is that the boundary conditions are different. For the cavity the boundary conditions give you a discrete density of states and standing waves while in free space you get a continuum of modes - some travelling left and some right.
$b_{in}(t)$ sort of represents the full solution to the Helmholtz equation in question while the $b(\omega)$ are the individual frequency components of $b_{in}(t)$.
The interesting bit is that in addition to asking about the total energy in the electromagnetic field (the Hamiltonian integral above) one can ask about the intensity or flux of energy passing through a surface. Now that we can talk about flux we can think of the $s^{-1}$ unit on $b^{\dagger}_{in}(t)b_{in}(t)$ as being related to the energy passing through a surface per unit time. 
For example, if we detect the output field $b_{out}(t)$ on a photodetector (which we often do) then we need only consider the energy in $b_{out}(t)$ in a narrow bandwidth (the detector response band). In that case the intensity is
$$
I(t) = \frac{c\epsilon_0}{2} |E(t)|^2 = \frac{c \epsilon_0}{2} \frac{\hbar \omega}{2V \epsilon_0} b_{out}^{\dagger}(t) b_{out}(t) = \frac{\hbar \omega c}{4 V} b_{out}^{\dagger}(t)b_{out}(t)
$$
We see that the units work out. There is the pesky mode volume $V$ showing up which came from quantizing the electric field in the continuum.. if things are handled properly one can make this go away. But the point is the factor $\hbar \omega$ gives us energy in the formula, the factor $\frac{c}{V}$ tells us something about the geometry of the travelling wave. $V$ might be related to the cross section of the beam. Then with all of this it is reasonable to say that $b_{out}^{\dagger}(t)b_{out}(t)$ is in fact telling us about the number of photons falling on the detector per unit time.
Anyways, this answer has sort of gotten away from me. My advice would be the following. First, realize that this formalism is working in the Heisenberg picture. Get comfortable with that. To help get comfortable with that it may be beneficial to try to translate everything into a purely classical description. In particular there is in fact very little quantum about the input output formalism except the fact that the operators all have funny commutation relations, so stripping away the quantumness might help.
Next I'll just give a list of all the things that might be confusing here having to do with the relationship between frequency and time.


*

*Heisenberg picture

*Fourier Transform

*Power spectral density

*amplitude vs energy/power

*energy vs intensity

*discrete vs continuous density of states


Anyways, I hope this is able to help a little bit. Sorry it's so scatterbrained. 
