Hermitian operator
Your definition is ok and equivalent to
$\langle a| A b \rangle = \langle A a| b \rangle$
Eigenspace of a Hermitian operator
It's possible to prove that a Hermitian operator has eigenvectors that form an orthogonal vector base, associated with real eigenvalues
$A | n \rangle = \lambda_n | n \rangle$.
Than you can normalize eigenvectors to get a orthonormal basis.
Generic wave function
Given a generic wave function $| \psi \rangle$ you can write it using the basis $\{ |n\rangle\}_n$, as a superposition of states
$| \psi \rangle = a_n |n\rangle $
so that you can write the action of the operator $A$ as
$ A | \psi \rangle = a_n A |n\rangle = a_n \lambda_n |n\rangle $.
Normalization condition and interpretation as a probability
With the normalization condition for the wave function $\langle \psi|\psi \rangle = 1$ (holding also for eigenvectors), you readily get
$1 =\langle \psi|\psi \rangle = \langle a_n n | a_m m \rangle = \sum_n a_n^2$.
If a system is ina state descirbe by wavefunction $|\psi \rangle$, the probability of measuring state $| m \rangle$ is
$|\langle \psi | m \rangle|^2 = |\langle a_n n | m \rangle|^2 = a_m^2$.
Eigenvalues and result of measurement process
On the other hand, eigenvalues of an operator are connected to the values you get from measurement process for the physical quantity associated with that oeprator. See Born's rule, https://en.m.wikipedia.org/wiki/Born_rule
There is a general argument for the existence of the unitary operator. From the Stone-von Neumann theorem, if you know where $a$ is mapped, then there is only one unitary operator that can represent the transformation (up to a global phase).
You can construct this operator more explicitly. One way is to realize that your transformations form a path connected Lie group. Formally, you can view it as the subgroup of matrices of size $2N\times2N$ of the form:
$$
\begin{pmatrix}
A & B\\
B^* & A^*
\end{pmatrix}
$$
with your additional constraints:
$$
A A^\dagger -BB^\dagger =1\\
AB^T-BA^T =0
$$
You can consider a path from the identity to your transformation, namely $A_t,B_t$ with $t\in[0,1]$ and:
$$
\begin{align}
A_0&=1 & B_0&=0 \\
A_1&=A& B_1&=B
\end{align}
$$
This gives a varying $a_t$ satisfying the ODE:
$$
\dot a_t= \alpha_ta_t+ \beta_ta_t^\dagger
$$
with the matrices $\alpha_t,\beta_t$ defined by:
$$
\frac{d}{dt} \begin{pmatrix}
A_t & B _t\\
B _t ^* & A _t ^*
\end{pmatrix} = \begin{pmatrix}
\alpha_t& \beta_t\\
\beta_t^* & \alpha_t^*
\end{pmatrix} \begin{pmatrix}
A_t & B _t\\
B _t ^* & A _t ^*
\end{pmatrix}
$$
You want to recognize the ODE as a Heisenberg equation from a suitably chosen quadratic Hamiltonian:
$$
\dot a_t =-i[a,H]
$$
Using the CCR’s you can check that a suitable candidate (unique up to additional multiple of identity) is:
$$
H_t= i(\alpha_t)_{ij}a_i^\dagger a_j+\frac{i}{2}(\beta_t)_{ij}a_i^\dagger a_j^\dagger-\frac{i}{2}(\beta_t^*)_{ij}a_i a_j
$$
$H_t$ is hermitian and gives the correct equations of motion using the relations (taking the derivatives of the constraints):
$$
\alpha_t+\alpha_t^\dagger =0 \\
\beta_t-\beta_t^T=0
$$
Your unitary operator is given by the time ordered exponential:
$$
U=\mathcal T \exp\left(\int_0^1dt H_t\right)
$$
Actually, it turns out that the exponential is surjective on the group, so you can choose
Choose a path with constant $\alpha_t,\beta_t$. This gives a constant $H_t$ and the time ordered exponential is a simple exponential.
Physically, you are constructing the unitary operator from harmonic oscillators and squeeze operators. The process is explicit in the case $N=1$.
Hope this helps.
Best Answer
My answers to questions 1-3 are as follows:
The argumentation is as follows.
Breuer and Petruccione define the the first standard form of the Liuvillian as \begin{equation} \mathcal{D}(\rho) = \sum_{ij} a_{ij} \left(F_i\rho F_j^\dagger - \frac{1}{2}\{\rho, F_j^\dagger F_i\}\right). \end{equation} It is not unique, as it is invariant under the following transformations of the operators $F_i$ and the matrix $a_{ij}$: \begin{equation} \begin{gathered} F_i = u_{ki} F'_k,\\ a_{kl}' = u_{ki} a_{ij}u^*_{lj} , \end{gathered} \end{equation} where $u_{ki}$ is an arbitrary unitary matrix.
The question is about the non--uniqueness of the diagonal form of the Liuvillian. The latter is obtained from the first standard form by choosing the unitary $u$ which diagonalizes the matrix $a_{ij}$: \begin{equation} u_{ki} u^*_{lj} a_{ij} = \Gamma_k^A\delta_{kl}, \end{equation} %which also transforms the set of operators $F_i$ into $A_i$: \begin{equation} F_i = u_{ki} A_k. \end{equation} The choice of the operators $A_k$ is also not unique: there always exist unitary matrices $U_{ij}$ which preserve the diagonal form of the Liuvillian after the transformation %and sets of $B_i$ \begin{equation} A_i = U_{ij}B_j. \end{equation} %which preserves the diagonal form of the Liuvillian and therefore the matrix $a_{ij} = \Gamma_i^A\delta_{ij}$: %Using the transformation law for the matrix $a_{ij}$, one gets the condition The condition for preserving the diagonal form follows from the transformation law for the matrix $a_{ij}$: \begin{equation} \sum_{k,l} U_{ik} \Gamma_k^A U_{j,k}^* = \Gamma_i^B\delta_{ij}. \end{equation} The class of unitaries $U$ which turns the diagonal matrix $\Gamma_{k}^A\delta_{kl}$ into a diagonal matrix is quite restricted. Such unitaries include:
Therefore, the sets of $\Gamma_k^A$ and $\Gamma_k^B$ differ only by permutation, whereas the operators $B_i$ from the question may be nontrivial linear combinations of $A_i$ providing that some of $\Gamma_k^A$ coincide.
Also, $U_{ij}$ is nonzero only providing $\Gamma_i^A = \Gamma_j^B$. Therefore, the transformation $A_i = U_{ij}B_j$ is equivalent to $\sqrt{\Gamma_i^A}A_i = U_{ij}\sqrt{\Gamma_j^B}B_j$ which was in the question.