Monopole operator: correlation functions Let's consider free Maxwell theory:
$$
S_{Maxwell} = \int d^dx \; -\frac{1}{4e^2}F_{\mu\nu}F^{\mu\nu}
$$
In such theory one can define monopole operator using path integral via correlation functions with other operators:
$$
\langle\mathcal{M}(x) \mathcal{O}_1(x_1)\dots \mathcal{O}_n(x_n)\rangle
=
\int D A_{\mu}\;\mathcal{O}_1(x_1)\dots \mathcal{O}_n(x_n) \;e^{-S_{Maxwell}}
$$
Here path integration goes through configurations, where in point x included magnetic flux across surface, that surrounds point x of insertion of monopole.
I wanna to understand, how calculate correlation functions between monopole operators (in arbitrary d!).
For example:
$$
\langle\mathcal{M}(x)\rangle = ???
$$
$$
\langle\mathcal{M}(x)\mathcal{M}(y)\rangle = ???
$$
How to calculate such correlators? I think that simplest example in Maxwell theory in 3d (monopole in such theory rather instantons), due to duality with compact scalar.
I will be very appreciated for any suggestions/answers!
Edit(3D case):
Following section 8.1.1 David Tong: Lectures on Gauge Theory:
$$
\langle\mathcal{M}(x)\rangle = \int D\sigma \;e^{i\sigma(x)}\; e^{-\int d^3x\; \frac{e^2}{8\pi^2}\partial_\mu \sigma \partial^\mu \sigma}
$$
$$
\langle\mathcal{M}(x)\mathcal{M}(y)\rangle =
\int D\sigma \;e^{i\sigma(x)}e^{i\sigma(y)} \;e^{-\int d^3x\; \frac{e^2}{8\pi^2}\partial_\mu \sigma \partial^\mu \sigma}
$$
Taking this Gaussian integrals using:
$$
\int D\sigma \;e^{i\int d^3x\;J\sigma}\; e^{-\int d^3x\; \frac{e^2}{8\pi^2}\partial_\mu \sigma \partial^\mu \sigma}
=
e^{\frac{2\pi^2}{e^2}\int d^3x d^3y\; J(x)\Box^{-1}(x-y)J(y)}
$$
I obtained:
$$
\langle\mathcal{M}(x)\rangle = e^{\frac{2\pi^2}{e^2}\Box^{-1}(0)}
$$
$$
\langle\mathcal{M}(x)\mathcal{M}(y)\rangle =e^{\frac{4\pi^2}{e^2}\Box^{-1}(0)} e^{\frac{4\pi^2}{e^2} \Box^{-1}(x-y)}
$$
$$
\Box^{-1}(x-y) = \frac{1}{|x-y|}
$$
Maxwell theory in 3d is scale invariant.. Why such correlators doesn't have such property?
Is my calculations true? How to generalize this to higher dimmensions?
 A: “Monopole operators” only exist as local operators in three dimensions. For free Maxwell theory in that dimension an easy way to think about them is in terms of the dual compact scalar, as you describe in your edit.
In four dimensions, instead of local monopole operators you have ‘t Hooft lines, and so on and so forth.
Edit: Regarding your calculations, you have to be careful. As explained in a different answer, the monopoles are charged under a $U(1)_T$ global symmetry ($T$ for topological). The monopole of charge $n$ is $e^{in\sigma}$.  The corresponding Ward identity implies that the only non-zero correlators are those for which the charges add up to zero. You calculations are in the right direction, but the symmetry forbids a nonzero one-point function (unless it is spontaneously broken). If you do the two-point function correctly you should see that the singular factors of $\Box^{-1}(0)$ cancel in the exponent due to charge conservation.
Also, Maxwell theory in $D=3$ is not scale invariant. The gauge coupling in $D\neq 4$ is always dimensionful.
A: As @jpm correctly says, you can have magnetic monopoles as local operators only in $3d$. The reason is that in $3d$ they are charged under a Noether current
$$
j^\mu = \tfrac{i}{2}\epsilon^{\mu\nu\rho} F_{\nu\rho}\,, \qquad \mbox{or} \qquad {}^\star j = F\,.
$$
Whenever you have a current, you also have a particle because the current naturally couples to its worldline
$$
S = \int \mathrm{d}\tau \,j\,.
$$
And in dimensions bigger than three of course you cannot write that because the $\epsilon$ does not have three indices.
But you can generalize the same concept to higher form symmetries $[1]$. Instead of having a $d-1$ form $\epsilon^{\mu_1\mu_2\ldots \mu_{d}} j_{\mu_1}$ you have a $d-p-1$ form ${}^\star j$. This form naturally couples to a $p+1$ dimensional worldvolume. So if $p=1$ it couples to the worldsheet of a line operator and so on.
These are called "symmetries" because you can define a topological operator by integrating them over a $d-p-1$ dimensional sphere
$$
Q(\Sigma) = \int_{\Sigma} {}^\star j\,.
$$
This operator plays the same role as a charge for ordinary symmetries. The catch is that the groups associated to higher form symmetries can only be abelian. For the case of monopoles indeed we have $U(1)$.

$\;[1]\;$ D. Gaiotto, A. Kapustin, N. Seiberg, B. Willett, Generalized Global Symmetries, ArXiv:1412.5148 [hep-th]
A: One way to define monopole operators, or their higher-dimensional generalisations,'t Hooft operators, as @jpm and @MannyC pointed out is to dualise the photon, $A$ to the magnetic photon, $\tilde{A}$ through
$$ \mathrm{d}\tilde{A} := \star \mathrm{d}A.\label{1}\tag{1}$$
Then define the 't Hooft operators as
$$ \text{'tH}(X) := \exp\!\bigg(i\int_X \tilde{A}\bigg).$$
These are like Wilson operators for the magnetic photon. Note that since the photon is a 1-form, the magnetic photon is a $(d-3)$-form, by counting degrees in (\ref{1}). So you have to integrate it on a $(d-3)$-manifold. Therefore $\text{'tH}(X)$ is extended in $(d-3)$-dimensions, confirming what both @jpm and @MannyC explain. It is a local operator in 3 dimensions, a line in 4 dimensions and so on.
Now if you want to calculate $n$-point functions of $\text{'tH}(X)$s you can dualise the whole lagrangian, so you have
$$\big\langle\text{'tH}(X_1)\cdots\text{'tH}(X_n)\big\rangle = \int \mathrm{D}\tilde{A}\ \text{'tH}(X_1)\cdots\text{'tH}(X_n) \exp\!\bigg(-\frac{1}{2e^2}\int \mathrm{d}\tilde{A}\wedge\star\mathrm{d}\tilde{A}\bigg)$$
and it looks like a Wilson line $n$-point function calculation, but for the magnetic photon.
However, if you want to calculate correlation functions involving both 't Hooft operators and local operators involving the normal photon, it is not that straightforward because these operators are usually not easily expressible in terms of the magnetic photon.
A: Let me try to answer your question.
As was rightfully mentioned in other responses 3d theories are the only theories when a monopole can be a local operator.
Thus, let us stick to a 3d case.
The Maxwell theory in 3d is dual to a free periodic scalar.
What does it mean?
Let's forget for a second that we have the Maxwell theory.
In 3d we do have a spontaneous symmetry breaking so our theory would start its life in a superselection sector labeled by $\phi$.
For a periodic scalar, the right operator is $M = e^{i \phi}$ as it's well defined.
So, let's choose $M$ to be equal to some number $e^{i a}$.
Notice that we can make it an arbitrary number on a circle $S^1$. Moreover, you can see that we don't actually calculate anything. We "choose" $M$ to be something knowing that there is an SSB and the IR structure of all possible vacua is recovered just from a classical picture.
That is an "extra" input for our path integral.
What comes next?
We can easily calculate
\begin{equation}
\left<M(x)M(y) \right> = e ^{ 2 i a + \frac{4\pi^2}{e^2} \frac{1}{|x-y|}}
\end{equation}
Going back to the Maxwell theory we see that the correlator stays the same with $M$ being a monopole operator.
