First of all, what's the motivation for the Yang-Mills action and how
should I understand the coupling constants $\theta$ and $g$?
I would say motivation comes from experiments. For instance it is an experimental fact that the electric charge is conserved. The associated current is also conserved, in the sense of
$$\partial_\mu J^\mu=0.$$
Therefore we can write this current as the curl of a vector potential $A^\mu$. Since the curl of the grad vanish, there is a redundancy
$$A^\mu(x)\rightarrow A^\mu(x)-\partial_\mu\Lambda(x),$$
giving the same measured current. This redundancy is called gauge symmetry and in the present case it is just $U(1)$. It happens that for some fundamental interactions there are more vector potentials, more charges and we end up with a gauge symmetry based on a non-Abelian group. This is a Yang-Mills theory. For example, the Quantum Chromodynamics, which describes the strong interaction and is based on the $SU(3)$ color symmetry.
The classical fields of the theory must be solutions of equations of motions which are obtained from the action, according to the Hamilton Principle. The quantum fields also require an action in order to quantize it using the path integral method. When constructing the action you have to show how the fields interact. The strength of these interactions are given by the magnitude of the ecoupling constants, which is an experimental input.
How can I get this so called expansion in power series with variable
$\tau$ of the probability amplitude?
Going from a non-interacting quantum field theory to an interacting one is in same sense similar to going from a harmonic oscillator to an anharmonic one. For example, you add a quartic term in the potential, for both cases. Then you solve the equation of motion of the oscillator or the probabilities amplitudes for the fields using perturbation theory. The probabilities in the quantum theory may be obtained from something called generating functional, which involves the exponential the action. Perturbative expansion here means to expand this exponential in powers either of $\hbar$ or of the coupling constants.
What was the motivation to start looking at this duality? A creation
of an everywhere defined (in $\tau$) gauge theory, maybe?
Some context: It was noted since the 70's that some non-Abelian gauge theories admit solutions with stable magnetic charge. Then Montonem and Olive and Goddard, Nuyts and Olive , noticed that they could map the mass spectrum of electric charges of a particular theory to the mass spectrum of magnetic charges of another theory (called dual), as long as they assume a particular map between the couplings of these theories. They conjectured that these theories are electromagnetic dual. This is a non-Abelian generalization of the electromagnetic duality in Maxwell theory. A particular example where this duality is conjectured is between two "copies" of Georgi-Glashow model: $SO(3)\rightarrow SO(2)$ with the Higgs in the adjoint.
$SL(2,\mathbb Z)$ duality: It was shown by Witten that the addition of a topological term ($\theta$-term) to the Georgi-Glashow model gives the spectrum of electric charges
$$q_e=e\left(n_e+\frac{\theta}{2\pi}n_m\right),\quad n_e,n_m\in\mathbb Z,$$
This theory also admit dyons, particles with the above electric charge and also a magnetic charge
$$q_m=\frac{4\pi}{e}n_m.$$
Then if we define
$$\tau=\frac{\theta}{2\pi}+\frac{4\pi i}{e^2},$$
the charge of the dyon can be written as
$$\mathcal Q = q_e+iq_m=e(n_e+\tau n_m).$$
The spectrum of the theory belongs then to a lattice whose vertices give the charges $(n_m,n_e)$. The masses of the dyons (of a theory with coupling $\tau$) are given by
$$M(n_e,n_m;\tau)=ve|n_e+n_m|=ve
\left|
\left(\begin{array}{cc}
n_m&n_e
\end{array}\right)
\left(\begin{array}{c}
\sqrt u\tau\\
\sqrt u
\end{array}\right)
\right|,$$
where $\sqrt u\equiv ev$ and $v$ is the Higgs vev. In order for a duality happen there must be a mapping between the mass spectrum of two theories, i.e.
\begin{equation}\label{sl6}
\left|
\left(\begin{array}{cc}
n_m&n_e
\end{array}\right)
\left(\begin{array}{c}
\sqrt u\tau\\
\sqrt u
\end{array}\right)
\right|
=
\left|
\left(\begin{array}{cc}
n'_m&n'_e
\end{array}\right)
\left(\begin{array}{c}
\sqrt{u'}\tau'\\
\sqrt{u'}
\end{array}\right)
\right|.
\end{equation}
A possible solution is
\begin{align*}
\left(\begin{array}{c}
\sqrt{u'}\tau'\\
\sqrt{u'}
\end{array}\right)
&=e^{i\varphi}\mathcal M
\left(\begin{array}{c}
\sqrt u\tau\\
\sqrt u
\end{array}\right),\\
\left(\begin{array}{cc}
n'_m&n'_e
\end{array}\right)
&=
\left(\begin{array}{cc}
n_m&n_e
\end{array}\right)\mathcal M^{-1},
\end{align*}
with $\varphi\in\mathbb R$ and
$$
\mathcal M=
\left(
\begin{array}{cc}
A&B\\
C&D
\end{array}
\right),
\quad \det \mathcal M= 1, \quad A,B,C,D\in\mathbb Z.
$$
Hence
$$\mathcal M\in SL(2,\mathbb Z).$$
The group $SL(2,\mathbb Z)$ has two generators
$$T=
\left(
\begin{array}{cc}
1&1\\
0&1
\end{array}
\right),
\quad
S=
\left(
\begin{array}{cr}
0&-1\\
1&0
\end{array}
\right).
$$
The duality transformations generated by $T$ and $S$ are called $T$-duality and $S$-duality, respectively. As you can see from the $S$ generator, the $S$-duality invert the coupling constant,
$$S:\tau\rightarrow-\frac 1\tau,$$
i.e. it is a duality that maps strong coupling to weak coupling. On the other hand $T$ acts as
$$T:\tau\rightarrow\tau+1,$$
which represents an invariance with respect of $\theta\rightarrow\theta+2\pi$.
It is important to notice that a possible electromagnetic duality is given only by a subgroup of $SL(2,\mathbb Z)$ because it is necessary to consider other quantum number of the particles.
The importance of the $S$-Duality is in the fact that you can use results obtained in a theory with weak coupling (where perturbation theory is valid) in the dual theory which has strong coupling (and perturbation theory is not valid).
