Question 1 is answered in parts 1 through to 6: the idea is that each part slowly builds the tools required to derive the process equation for $S_t$ under the $S_t$ Numeraire.
Question 2 & Question 3 are then answered in part 7.
- Part 1: Expectation of a function of a Random variable:
Let $X(t)$ be some generic Random Variable with probability density function given by $f_{X_t}(h)$, where $h$ is a "dummy" variable. Let $g(X_t)$ be some (well-behaved) function of $X_t$. Then (I am stating the below without proof):
$$\mathbb{E}[g(X_t)]=\int_{-\infty}^{\infty}g(X_t)f_{X_t}(h)dh$$
- Part 2: Radon-Nikodym Derivative:
Let $\mathbb{P^1}$ be a Probability measure defined via the Probability Density Function of some random variable $X_t$:
$$\mathbb{P^1}(A):=\int_{-\infty}^{a}f_{X_t}(h)dh$$
For all events $\{A: X_t \leq a\}$.
Radon-Nikodym derivative is implicitly defined as some Random-Variable (let's call it $Y_t$) that satisfies the following:
$$ \mathbb{P^2}(A) = \mathbb{E^{P^1}}[Y_t \mathbb{I_{\{ A\}}}] $$.
The above definition becomes more intuitive with a specific example: let $X_t$ be a standard Brownian Motion, i.e. $X_t:=W_t$, and let $Y_t:=e^{-0.5\sigma^2t+\sigma W_t}$. Basically $Y_t=g(W_t)$, where $g()$ is a well-behaved function: so we can make use of the result in part 1, specifically:
$$ \mathbb{E^{P^1}}[Y_t \mathbb{I_{\{ A\}}}] = \mathbb{E^{P^1}}[g(W_t) \mathbb{I_{\{ A\}}}] = \\ = \int_{-\infty}^{\infty}g(X_t)f_{X_t}(h) \mathbb{I_{ \{ W_t \leq a \}}}dh = \\ = \int_{-\infty}^{a}g(X_t)f_{X_t}(h)dh = \\ = \int_{-\infty}^{a}e^{-0.5\sigma^2t+\sigma h}\frac{1}{\sqrt{2\pi}}e^{\frac{-h^2}{2t}}dh = \\ =\int_{h=-\infty}^{h=k}\frac{1}{\sqrt{2\pi}}e^{\frac{-(h^2-\sigma t)}{2t}}dh $$
(To go from the penultimate line to the last line, we just need to complete the square).
The main point: by applying the definition $\mathbb{P^2}(A) = \mathbb{E^{P^1}}[Y_t \mathbb{I_{\{ A\}}}]$, we can see how $Y_t$ "creates" a new probability measure: under $\mathbb{P^2}$, the same event, specifically $A: W_t \leq a$ has an altered probability, compared to the same event under $\mathbb{P^1}$.
By inspecting the probability $\mathbb{P^2}(A)=\mathbb{P^2}(W_t \leq a) = \int_{h=-\infty}^{h=k}\frac{1}{\sqrt{2\pi}}e^{\frac{-(h^2-\sigma t)}{2t}}dh$, we ca see that what was standard Brownian motion under $\mathbb{P^1}$ now has a probability distribution of a Brownian motion with a drift: so under $\mathbb{P^2}$, $W_t$ is no longer a standard Brownian motion, but a Brownian motion with drift $\sigma t$.
- Part 3: Cameron-Martin-Girsanov Theorem:
The theorem states that:
If $W_t$ is standard Brownian motion under some $\mathbb{P^1}$, then there exists some $\mathbb{P^2}$ under which $W_t$ is a Brownian motion with drift $\mu t$. The Radon-Nikodym derivative to get us from $\mathbb{P^1}$ to $\mathbb{P^2}$ is:
$$ \frac{d \mathbb{P^2}}{d \mathbb{P^1}}(t)= e^{-0.5\mu^2t+\mu W_t}$$
If $\tilde{W_t}:=W_t + \mu t$ is a Brownian motion with some drift $\mu t$ under some $\mathbb{P^1}$, then there exists some $\mathbb{P^2}$ under which $\tilde{W_t}$ is a standard Brownian motion (i.e. no drift). The Radon-Nikodym derivative to get us from $\mathbb{P^1}$ to $\mathbb{P^2}$ is:
$$ \frac{d \mathbb{P^2}}{d \mathbb{P^1}}(t)= e^{+0.5\mu^2t-\mu W_t}$$
We basically "proved" the C-M-G theorem in part 2 above.
- Part 4: Numeraire and Probability Measures
Under the risk-neutral measure, with deterministic money market as Numeraire, the stock price process is: $S_t=S_0exp\left[ (r-0.5 \sigma^2)t+\sigma W(t) \right]$. The only source of randomness in this process is $W_t$, which is a standard Brownian motion under $\mathbb{P^Q}$ associated with the Numeraire $N_t:=e^{rt}$.
Since $W_t$ is the only source of randomness, this gives us an idea of how a change of probability measure will work for the process $S_t$: the change of measure will be driven via a Radon-Nikodym derivative applied to $W_t$. If we can somehow get a Radon-Nikodym derivative that resembles the one from the C-M-G Theorem, then we're in for an easy change of measure: we could apply the CMG theorem directly to $W_t$ in the process equation for $S_t$!!
- Part 5: Change of Numeraire formula
Without proof, if we want to change numeraire from $N_t$ to some $N^{2}_t$, the Radon-Nikodym derivative we need to use is:
$$ \frac{dN^{2}_t}{dN_t}:= \frac{N(t_0)N_2(t)}{N(t)N_2(t_0)} $$
(The proof of the above formula can be found here: Change of Numeraire formula)
- Part 6: Choosing $S_t$ as Numeraire
Applying the formula from part 5 above, we get:
$$ \frac{dN^{S_t}_t}{dN_t}:= \frac{N(t_0)N^{S_t}(t)}{N(t)N^{S_t}(t_0)} = \\= \frac{1*S_t}{e^{rt}S_0}= \\ = \frac{S_0\exp\left[ (r-0.5 \sigma^2)t+\sigma W(t) \right]}{e^{rt}S_0}= e^{-0.5\sigma^2t+\sigma W_t} $$
The above result is great news, because we can use part 3 directly and apply $e^{-0.5\sigma^2t+\sigma W_t}$ as Radon-Nikodym derivative to $W_t$: we know this will introduce the drift $\sigma t$ under the probability measure defined through $\frac{dN^{S_t}_t}{dN_t}=e^{-0.5\sigma^2t+\sigma W_t}$.
Let $\tilde{W_t}:=W_t-\sigma t$ be a Brownian motion with a drift equal to $-\sigma t$ under $\mathbb{P^Q}$. Inserting $\tilde{W_t}$ into the process equation for $S_t$ under $\mathbb{P^Q}$, we get (pure algebraic manipulation, no tricks here):
$$S_t=S_0\exp\left[ (r-0.5 \sigma^2)t+\sigma W(t) \right]= \\ = S_0\exp\left[ (r-0.5 \sigma^2)t+\sigma (\tilde{W}(t)+\sigma t) \right] = \\ = S_0\exp\left[ (r-0.5 \sigma^2)t+\sigma^2 t + \tilde{W}(t) \right] = \\ = S_0\exp\left[ (r+0.5 \sigma^2)t+ \tilde{W}(t) \right]$$
The above equation is not particularly useful in any way. But we can now do the following: we can apply the Cameron-Martin-Girsanov theorem to $\tilde{W}_t$, which is very convenient: taking the Radon-Nikodym drivative $\frac{dN^{S_t}_t}{dN_t}=e^{-0.5\sigma^2t+\sigma W_t}$ and applying it to $\tilde{W_t}$ will add the drift $\sigma t$. But $\tilde{W_t}$ has negative drift equal to $-\sigma t$. Therefore, the Radon-Nikodym derivative $\frac{dN^{S_t}_t}{dN_t}$ will "kill" the drift of $\tilde{W_t}$. Consequently, under the probability measure associated with $S_t$ as Numeraire, $\tilde{W_t}$ becomes a standard Brownian motion with no drift.
That's why under the Stock numeraire, the process for the stock price becomes (with $\tilde{W}_t$ being a standard Brownian motion):
$$S_t=S_0\exp\left[ (r+0.5 \sigma^2)t+ \tilde{W}(t) \right]$$
It is worth noting that often people use "lazy" notation and don't put the 'tilde' sign on the Brownian motion under the new measure: but I prefer to do it to emphasize that it's a different process to the plain Brownian motion $W_t$ under the risk-neutral measure.
Part 7: evaluating $\mathbb{E^{N_{S}}}[S_t\mathbb{I_{\{S_t > k\}}}]$:
I think there are multiple ways the expectation can be evaluated. The method that uses the least advanced mathematics but involves the most labor is direct evaluation via an integral:
$$ \mathbb{E^{N_{S}}}[S_t\mathbb{I_{\{S_t > k\}}}] = \int_{S_t=k}^{\infty} S_t f_{S_t}(S_t)dS_t = \int_{h=k}^{\infty} h f_{S_t}(h)dh $$
We know that $S_t$ is log-normally distributed, so we know the density of $S_t$ (https://en.wikipedia.org/wiki/Log-normal_distribution):
$$f_{S_t}(h)= \frac{1} {h \sqrt{t}\sigma \sqrt{2\pi}} e^{-\frac{(ln(h/S_0)-(r-0.5\sigma^2)t)^2}{2\sigma^2t}}$$
Plugging this into the integral results in the cancellation of the $h$ in the first denominator:
$$\int_{h=k}^{\infty} \frac{1} {\sqrt{t}\sigma \sqrt{2\pi}} e^{-\frac{(ln(h/S_0)-(r-0.5\sigma^2)t)^2}{2\sigma^2t}}dh $$
I am gonna do the following substitutions: $y:=ln(h/S_0)$, so that $h=S_0e^e$, $dh=S_0e^ydy$, and when $h=K$, we get $y=ln\left( \frac{K}{S_0} \right)$.
Integrating via substitution then yields:
$$\int_{y=ln(K/S_0)}^{\infty} \frac{1}{\sigma \sqrt{t}} \frac{1}{\sqrt{2 \pi}} e^{\frac{(y-(r-0.5\sigma^2)t)^2}{2\sigma^2t}}S_0 e^y dy$$
I am now gonna simplify the notation further with: $\tilde{\mu}:=(r-0.5\sigma^2)t$ and $\tilde{\sigma}:=\sigma \sqrt{t}$, so the integral becomes:
$$\int_{y=ln(K/S_0)}^{\infty} \frac{1}{\tilde{\sigma}} \frac{1}{\sqrt{2 \pi}} e^{\frac{(y-\tilde{\mu})^2}{2\tilde{\sigma}^2}}S_0 e^y dy$$
Completing the square between $e^y$ and $e^{\frac{(y-\tilde{\mu})^2}{2\tilde{\sigma}^2}}$ gives:
$$ \exp(y) \exp\left(\frac{(y-\tilde{\mu})^2}{2\tilde{\sigma}^2}\right) = \\ = \exp \left(\frac{(y-(\tilde{\mu}+\tilde{\sigma}))^2}{2\tilde{\sigma}^2}\right)*\exp\left(\tilde{\mu}+0.5\tilde{\sigma}^2\right) = \\ =\exp \left(\frac{(y-(\tilde{\mu}+\tilde{\sigma}))^2}{2\tilde{\sigma}^2}\right)*\exp\left(rt\right) $$
The last line uses the fact that $\tilde{\mu}+0.5\tilde{\sigma}^2=(rt-0.5\sigma^2t)+0.5\sigma^2t=rt$.
Plugging back into the integral gives:
$$S_0e^{rt}\int_{y=ln(K/S_0)}^{\infty} \frac{1}{\tilde{\sigma}} \frac{1}{\sqrt{2 \pi}} \exp \left(\frac{(y-(\tilde{\mu}+\tilde{\sigma}))^2}{2\tilde{\sigma}^2}\right)dy$$
Finally, one last substitution: I will take $z:=\frac{y-(\tilde{\mu}+\tilde{\sigma}^2)}{\sqrt{t}\sigma}$, which gives $dy=\sqrt{t}\sigma dz$. Furthermore, when $y=ln\left( \frac{K}{S_0} \right)$, we get:
$$z=\frac{ln\left( \frac{K}{S_0} \right)-(\tilde{\mu}+\tilde{\sigma}^2)}{\sqrt{t}\sigma}=\frac{ln\left( \frac{K}{S_0} \right)-(rt+0.5 \sigma^2t)}{\sqrt{t}\sigma} = \\ = (-1) \frac{ln\left( \frac{S_0}{K} \right)+rt+0.5 \sigma^2t}{\sqrt{t}\sigma} = -d_1 $$
So plugging this last substitution for $y$ into the integral gives:
$$S_0e^{rt}\int_{y=ln(K/S_0)}^{\infty} \frac{1}{\tilde{\sigma}} \frac{1}{\sqrt{2 \pi}} exp \left(\frac{(y-(\tilde{\mu}+\tilde{\sigma}))^2}{2\tilde{\sigma}^2}\right)dy= \\ = S_0e^{rt}\int_{z=-d_1}^{\infty} \frac{1}{\sqrt{2 \pi}} exp \left(\frac{z^2}{2} \right)dz= \\ =S_0e^{rt}\mathbb{P}(Z>-d_1)=S_0e^{rt}\mathbb{P}(Z \leq d_1) = S_0e^{rt} N(d_1) $$